**Algorithms and Methods for Designing and Scheduling Smart Manufacturing Systems**

Editors

**Vladimir Modrak Zuzana Soltysova**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Vladimir Modrak Technical University of Kosice ˇ Slovakia

Zuzana Soltysova Technical University of Kosice ˇ Slovakia

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: https://www.mdpi.com/journal/applsci/special issues/Algorithms Methods Smart Manufacturing Systems).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-4509-7 (Hbk) ISBN 978-3-0365-4510-3 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


### **Gedas Baranauskas and Agota Giedr˙e Raiˇsien˙e**


## **About the Editors**

#### **Vladimir Modrak**

Prof. Dr. Vladimir Modrak is Full Professor of Manufacturing Technology at Faculty of Manufacturing Technologies. His research interests include manufacturing logistics, cellular manufacturing, lean manufacturing, mass customization and other related disciplines. He was the leading editor and co-author of the books "Operations Management Research and Cellular Manufacturing Systems: Innovative Methods and Approaches" (2012), "Handbook of Research on Design and Management of Lean Production Systems" (2014), "Mass Customized Manufacturing: Theoretical Concepts and Practical Approaches" (2016), "Industry 4.0 for SMEs Challenges, Opportunities and Requirements" (2020), and "Implementing Industry 4.0 in SMEs: Concepts, Examples and Applications" (2021). He is Fellow of the European Academy of Industrial Management (AIM).

#### **Zuzana Soltysova**

Dr. Zuzana Soltysova completed her bachelor, master and PhD. study at Technical University of Kosice, Faculty of Manufacturing Technologies. Currently, she is professor assistant at Technical University of Kosice, Faculty of Manufacturing Technologies, Department of Industrial Engineering and Informatics. Her doctoral thesis was titled "Research of product and production complexity in terms of mass customization". Moreover, her research activities include complexity problems, axiomatic design, product and process modularity and related disciplines to Manufacturing management.

### *Editorial* **Algorithms and Methods for Designing and Scheduling Smart Manufacturing Systems**

**Vladimir Modrak \* and Zuzana Soltysova \***

Faculty of Manufacturing Technologies, Technical University of Kosice, Bayerova 1, 080 01 Presov, Slovakia **\*** Correspondence: vladimir.modrak@tuke.sk (V.M.); zuzana.soltysova@tuke.sk (Z.S.)

#### **1. Introduction**

This Special Issue is a collection of some of the latest advancements in designing and scheduling smart manufacturing systems. The smart manufacturing concept is undoubtedly considered a paradigm shift in manufacturing technology. This conception is part of the Industry 4.0 strategy, or equivalent national policies, and brings new challenges and opportunities for the companies that are facing tough global competition. Industry 4.0 should not only be perceived as one of many possible strategies for manufacturing companies, but also as an important practice within organizations. The main focus of Industry 4.0 implementation is to combine production, information technology, and the internet [1]. Therefore, an introduction of smart manufacturing systems is primarily associated with the adaptation of the Internet of Things, cyber physical systems, artificial intelligence, advanced robotics, cloud technology, and so forth. In particular, web technologies act as enablers of smart manufacturing that can promote a disruptive innovation in small and medium enterprises. However, some recent studies (see, e.g., [2–4]) have shown that preexisting managerial methods and philosophies such as lean manufacturing, reconfigurable manufacturing systems, or cellular manufacturing systems are of utter importance for the concept of smart manufacturing. In this context, manufacturing system design and scheduling methods also play a vital role in the era of the fourth industrial revolution. It is particularly evident from the link between Industry 4.0 and lean manufacturing that both domains are strongly interrelated [5]. Thus, we hope that this Special Issue may be helpful in solving the particular problems which are related to smart manufacturing systems and advanced manufacturing technologies.

#### **2. Description of the Papers**

The presented Special Issue consists of ten research papers presenting the latest works in the field. The papers include various topics, which can be divided into three categories— (i) designing and scheduling manufacturing systems (seven articles), (ii) machining process optimization (two articles), (iii) digital insurance platforms (one article). Most of the mentioned research problems are solved in these articles by using genetic algorithms, the harmony search algorithm, the hybrid bat algorithm, the combined whale optimization algorithm, and other optimization and decision-making methods.

The above-mentioned groups of articles are briefly described in this order in the rest of this editorial paper.

The paper written by J.S. Park, H.Y. Ng, T.J. Chua, Y.T. Ng, and J.W. Kim [6] presents a novel genetic algorithm approach that utilizes a multiple chromosome scheme to solve the flexible job-shop scheduling problem which involves two kinds of decisions—machine selection and operation sequencing. This novel genetic algorithm approach enables the application of an identical crossover strategy in the categorical and sequential parts. The authors used this unified approach for the extension of the existing candidate order-based genetic algorithm to the unified candidate order-based genetic algorithm to solve flexible job-shop scheduling problems.

**Citation:** Modrak, V.; Soltysova, Z. Algorithms and Methods for Designing and Scheduling Smart Manufacturing Systems. *Appl. Sci.* **2022**, *12*, 3011. https://doi.org/ 10.3390/app12063011

Received: 10 March 2022 Accepted: 14 March 2022 Published: 16 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The other paper titled "Calibration of GA parameters for layout design optimization problems using design of experiments" [7] focuses on finding out how the solutions of the cell formation problem are influenced by a set of probability parameters of genetic operators, namely crossover and mutation, including balanced weight factors. In view of this, the presented work attempts to employ the Taguchi approach to find an optimal combination of parameters that impact the efficiency of the genetic algorithm and to explore whether the optimal combination of the genetic operators for the given type of machine-cell formation problems can be influenced by the magnitude of the noise factors.

The next paper [8] introduces a robust cluster algorithm based on the concept of group technology. The paper is written by L.M. Li and its originality lies in its use of the algorithm to balance an assembly line by matching operators to workstations so that the line's workstations achieve the same targeted output rates. It enables managers to solve the problem of worker absences by assigning more than one operator with the required skillset to each workstation and rearranging them as needed. This algorithm has been applied to cellular manufacturing system problems. Four examples were presented to implement and validate this algorithm, where training time was reduced by matching operators' training and skills according to the workstations' skill requirements.

The paper written by L. Nagarajan, S.K. Mahalingam, S. Salunkhe, E.A. Nasr, J.P. Davim, and H. Hussein [9] proposes a novel methodology for simultaneous minimization of manufacturing objectives in tolerance allocation of complex assembly tasks. The methodology consists of a two-step process. For this purpose, a heuristic approach was applied to determine the best machine for each process. Subsequently, it applies a combined whale optimization algorithm with a univariate search method to allocate optimum tolerances with the best process selection for each sub-stage/operation. The proposed methodology was validated by solving two typical tolerance allocation problems of complex assemblies: a wheel mounting assembly and a knuckle joint assembly. Authors also compared their approach with existing ones. The comparison results showed that the proposed approach considerably reduces tolerance cost and machining time.

The next paper which is authored by S.K Mahalingam, L. Nagarajan, S. Salunkhe, E.A. Nasr, J.P. Davim, and H. Hussein aims to acquire the maximum number of nonlinear assemblies with closer assembly tolerance specifications by mating the different bins' components [10]. Their proposed approach is based on using the combination of the univariate search method and the harmony search algorithm. The efficacy of the proposed method is demonstrated by showing 24.9% of cost savings while making overrunning clutch assembly in comparison with the existing method. Based on the obtained results in this work, the contribution of the proposed novel methodology is legitimate in solving selective assembly problems.

The originality of the following paper written by Zheng J. and Wang Y. [11] lies in an application of the hybrid bat optimization algorithm, which is based on the variable neighborhood structure, and two learning strategies to solve a three-stage distributed assembly permutation flow shop scheduling problem (DAPFSP) with the aim to minimize the makespan. The proposed algorithm is firstly designed to increase the population diversity by classifying the populations which solves the difficult trade-off between convergence and diversity of the bat algorithm. Secondly, a selection mechanism is used to update the bat's velocity and location, solving the difficulty of the algorithm trade-off of exploration and mining capacity. For this purpose, the Gaussian learning strategy and elite learning strategy assist the whole population in jumping out of the local optimal frontier. Based on the obtained simulation results, this algorithm can solve the DAPFSP problem well. In comparison with other metaheuristic algorithms, it provides better performance than the compared ones; thus, it is suitable to find the optimal solution(s).

The brief overview of the further paper titled "A multi-criteria assessment of manufacturing cell performance using the AHP method" [12] is as follows. It introduces the solution for how to find the optimal manufacturing cell design from alternative designs by using a multi-criteria evaluation. Alternative design solutions are mutually compared

by using the selected performance criteria, namely operational complexity, production line balancing rate, and makespan. Then, multi-criteria decision analysis based on the analytic hierarchy process method is used to show that two more-cell solutions better satisfy the determined criteria of manufacturing cell design performance than three less-cell solutions. The main benefit of this approach lies in the exactly enumerated values of the selected criteria, according to which the points from the mentioned scale are assigned to the alternatives.

The first paper of the second group of articles is titled "Meta-heuristic technique-based parametric optimization for electrochemical machining of Monel 400 alloys to investigate the material removal rate and the sludge" [13]. The authors investigate in this work the predominant electrochemical machining process parameters, namely applied voltage, flow rate, and electrolyte concentration, and their effects on the performance measures, i.e., the material removal rate and the nickel presence. For this purpose, authors used a metaheuristic algorithm named as the grey wolf optimizer for the multi-objective optimization of the process parameters for electrochemical machining. The obtained results were compared with the moth–flame optimization and particle swarm optimization algorithms. Based on the obtained results, it was observed that all the process variables significantly influenced the objectives. Then, it is confirmed that these metaheuristic algorithms—the moth–flame optimization algorithm and the grey wolf optimizer—are suitable for finding the optimum process parameters for machining Monel 400 alloys with electrochemical machining.

The second paper of this group is named "Optimization of process parameters for turning Hastelloy X under different machining environments using evolutionary algorithms: A comparative study" [14]. In their research, the authors investigated the machinability of turning Hastelloy X with a PVD Ti-Al-N coated insert tool in dry, wet, and cryogenic machining environments. The machinability indices, namely cutting force, surface roughness, and cutting temperature, are studied for the different set of input process parameters such as cutting speed, feed rate, and machining environment. They used the experiments conducted as per L27 orthogonal array. The authors proposed the moth–flame optimization to identify the optimal set of turning parameters through the multiple linear regression models, in view of minimizing the machinability indices. The effectiveness of the proposed algorithm is evaluated in comparison to the findings of genetic, Grass–Hooper, grey wolf, and particle swarm optimization algorithms. Based on the obtained results, this algorithm outperformed the others.

The last article in the given order is titled "Reflections on the customer decision-making process in the digital insurance platforms: An empirical study of the Baltic market" [15]. It aims to expand upon the existing scientific knowledge of end-user behavioral patterns and process frameworks in the Baltic insurance market, by including and examining a factor group of technological enablers. This paper is focused on research results in digitalization, personalization, and customization levels within the Baltic non-life insurance market in Estonia. There are also three major factor groups and process stages identified which influence insurance purchase decision-making in digital insurance platforms in the Baltic market.

#### **3. Conclusions**

It is believed that the collection of the ten papers in this Special Issue will be beneficial to readers who are interested in applying modern algorithms and methods for designing and scheduling smart manufacturing systems and related problems. Although this Special Issue has been closed, the current market challenges justify the need for an additional in-depth research in this domain.

**Author Contributions:** Conceptualization, V.M. and Z.S.; formal analysis, V.M. and Z.S.; writing original draft preparation, V.M. and Z.S.; writing—review and editing, V.M. and Z.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the KEGA Project No. 025TUKE-4/2020 granted by the Ministry of Education of the Slovak Republic.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** Moreover, we would like to thank all the authors who participated in this Special Issue and all the reviewers, and give many thanks to the editorial team of the Applied Sciences journal.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Unified Genetic Algorithm Approach for Solving Flexible Job-Shop Scheduling Problem**

**Jin-Sung Park 1, Huey-Yuen Ng 2, Tay-Jin Chua 2, Yen-Ting Ng <sup>2</sup> and Jun-Woo Kim 1,\***


**Abstract:** This paper proposes a novel genetic algorithm (GA) approach that utilizes a multichromosome to solve the flexible job-shop scheduling problem (FJSP), which involves two kinds of decisions: machine selection and operation sequencing. Typically, the former is represented by a string of categorical values, whereas the latter forms a sequence of operations. Consequently, the chromosome of conventional GAs for solving FJSP consists of a categorical part and a sequential part. Since these two parts are different from each other, different kinds of genetic operators are required to solve the FJSP using conventional GAs. In contrast, this paper proposes a unified GA approach that enables the application of an identical crossover strategy in both the categorical and sequential parts. In order to implement the unified approach, the sequential part is evolved by applying a candidate order-based GA (COGA), which can use traditional crossover strategies such as one-point or two-point crossovers. Such crossover strategies can also be used to evolve the categorical part. Thus, we can handle the categorical and sequential parts in an identical manner if identical crossover points are used for both. In this study, the unified approach was used to extend the existing COGA to a unified COGA (u-COGA), which can be used to solve FJSPs. Numerical experiments reveal that the u-COGA is useful for solving FJSPs with complex structures.

**Keywords:** flexible job-shop scheduling problem; combinatorial optimization; genetic algorithm; candidate order-based genetic algorithm; multichromosome

#### **1. Introduction**

Production scheduling is one of the most important decision-making procedures on manufacturing shopfloors, as it helps to utilize resources efficiently and maintain competitiveness in manufacturing companies [1–3]. Most production scheduling problems are NP-hard combinatorial optimization problems such that the optimal schedules are hard to find in polynomial time [4]. In this context, metaheuristic algorithms that can be used to find near optimal schedules in practical time are popular in production scheduling literature [5]. Examples of metaheuristic algorithms for solving production scheduling problems are the genetic algorithm (GA) [6], tabu search (TS) [7], simulated annealing (SA) [8], and particle swarm optimization (PSO) [9].

This paper aims to develop a novel GA for solving the flexible job-shop scheduling problem (FJSP). The classical job-shop scheduling problem (JSP) is one of the most wellknown production scheduling problems. It consists of *m* machines and *n* jobs. A job contains a number of operations to be processed in a fixed order. In the JSP, each operation can be processed by one specific machine [10,11]. The FJSP is an extension of the classical JSP, which allows an operation to be processed by one of two or more machines. In other words, an operation is processed by one of alternative machines in the FJSP [12]. Since a machine for a specific operation is predetermined, the JSP can be solved by specifying

**Citation:** Park, J.-S.; Ng, H.-Y.; Chua, T.-J.; Ng, Y.-T.; Kim, J.-W. Unified Genetic Algorithm Approach for Solving Flexible Job-Shop Scheduling Problem. *Appl. Sci.* **2021**, *11*, 6454. https://doi.org/10.3390/app11146454

Academic Editors: Paolo Renna, Vladimir Modrak and Zuzana Soltysova

Received: 31 May 2021 Accepted: 9 July 2021 Published: 13 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the priorities of given operations such that a high-priority operation precedes others with lower priorities in the queue for given machines. Thus, the JSP can be regarded as a sequencing problem [13,14]. In comparison with the JSP, the FJSP requires an additional decision related to which machine is used to process a specific operation. Consequently, two kinds of decisions, i.e., operation sequencing and machine selection, are required to solve the FJSP [1].

GA has been a popular metaheuristic algorithm for solving the JSP and FJSP in recent decades [10–12,15]. In order to apply GA, a solution of the given problem must be encoded in the form of a string, known as a chromosome [6]. One issue is that a chromosome for the FJSP must consist of two subchromosomes to account for the operation sequence (OS) part and the machine selection (MS) part. In other words, existing GAs for the FJSP are based on a multichromosome. Typically, the MS part is represented by a string of categorical values, whereas the OS part forms a sequence of given operations. Since the MS and OS parts have quite different structures, it is difficult to apply an identical genetic operator to both of them. Consequently, the MS and OS parts are separately evolved by different kinds of genetic operators in existing GAs for the FJSP. This traditional approach has the following limitations: First, it makes the GA hard to implement since different kinds of genetic operators should be applied. In particular, crossover operators for sequencing problems are a lot more complicated than simple one-point or two-point crossovers [16]. Second, there is no correlation between the order and machine in a single operation. In other words, the evolution of the MS part occurs independently from that of the OS part. Assume that an operation has the optimal order and is assigned to the optimal machine in a chromosome. Then, it is desirable that this combination of order and machine should also be maintained in the chromosome's offspring, which is hard to achieve using the traditional approach. Consequently, the traditional approach can cause a loss of good combinations in chromosomes and unnecessary diversity in the population. To this end, this paper poses two research questions: how can we reduce the number of genetic operators in the GA for the FJSP, which will have a significant impact on the complexity of the GA? Can the GA with a reduced number of genetic operators, designed to maintain good combinations of OS and MS genes, deal with the FJSP effectively?

In order to overcome the limitations of the traditional approach, this paper proposes a novel unified GA approach for the FJSP, which enables the application of an identical crossover operator to both the OS and MS parts. The main idea of the proposed approach is to apply a candidate order-based GA (COGA) based on an identical crossover point to both the OS and MS parts. The COGA is a type of GA developed for solving sequencing problems. The most distinguished feature of the COGA is that simple point-based crossover strategies, including one-point and two-point and uniform crossovers, can be implemented [14]. This enables the application of an identical point-based crossover strategy with an identical crossover point to the OS and MS parts, which may help to maintain the combinations of good order and good machine in parent solutions.

The remainder of this paper is organized as follows: In Section 2, a literature review concerning the GA for the FJSP and COGA is provided. Section 3 outlines the concept of the unified GA approach, and procedures and genetic operators of the unified COGA (u-COGA) are introduced. Section 4 presents the experimental results obtained by applying the u-COGA to various benchmark FJSPs. Finally, Section 5 concludes the paper.

#### **2. Research Background**

#### *2.1. Genetic Algorithm for Solving the Flexible Job-Shop Problem*

FJSP is a well-known extension of the JSP, a classical production scheduling problem. During the past decades, GAs have been widely used to solve production scheduling problems, such as the JSP and FJSP [15]. Typically, a GA chromosome for the JSP contains only OS-type genes, since a solution for the JSP can be generated using a sequence of given operations [10,11,13,17]. In contrast, a GA chromosome for the FJSP typically consists of two types of genes: MS genes and OS genes, which are manipulated by applying different

kinds of genetic operators. Genetic operators adopted by previous research papers on the GA for the FJSP are summarized in Table 1.

Gao et al. [18] used extended order crossover (OX) and uniform crossover for the OS and MS parts, respectively. Furthermore, the authors proposed two kinds of mutation operators: random alternative for the MS part and immigration for the OS part.

In Pezzella et al. [19], the MS part is manipulated by assignment crossover and two mutation operators: random alternative mutation and greedy mutation. Assignment crossover is designed to exchange machine assignment information for certain operations. In other words, only some of the MS genes in a chromosome participate in a single assignment crossover. Random alternative mutation is used to assign an operation to a machine randomly chosen from the set of all alternative machines that process the operation. On the contrary, greedy mutation assigns an operation assigned to the machine with the maximum workload to the machine with the minimum workload. On the contrary, precedence preserving order-based crossover (POX) and precedence preserving shift (PPS) mutation were applied to the OS part. Pazzella et al. [19] used a GA with the aforementioned operators to minimize the makespan of FJSP. Moreover, Defersha and Chen [20] utilized a similar GA for the purpose of minimization of the makespan of the FJSP with sequence-dependent setup times.

Lei [21] applied the GA to solve the FJSP with fuzzy processing times. A two-point crossover operator was adopted for the MS part, whereas two kinds of sequencing crossover operators, i.e., generalized position crossover (GPX) and generalization of the precedence preservative crossover (GPPX), were used to recombine the genes in the OS part. Swap mutation was used to randomly modify the OS part; however, the GA in Lei [21] did not contain the mutation procedure for the MS part.

Wang et al. [22] applied multipoint preservative crossover (MPX) to the MS part and improved precedence operation crossover (IPOX) to the OS part. They used random alternative mutation and greedy mutation to randomly modify the MS part in the chromosome. The authors proposed a greedy mutation operator designed to assign an operation to a machine that provides the shortest processing time (SPT). For the OS part, conventional insertion mutation was adopted.

Zhang et al. [1] used classical two-point crossover and uniform crossover operators to recombine the genes in the MS part of parent solutions, while the genes in the OS part are recombined using preserving order-based crossover (POX). As in Wang et al. [22], the genes in the MS part are mutated by greedy mutation based on an SPT machine. In contrast, the mutation for the OS part is performed by applying swap mutation.

Jiang and Du [23] used two-point crossover for the MS part and POX for the OS part. The mutation operations for the MS part and OS part were random alternative and swap mutation, respectively.

Türkyılmaz and Bulkan [24] adopted three types of crossover operators for solving the FJSP. Uniform crossover was used to perform the crossover operation for the MS part, while two variations of POX, modified-1 POX (MPX1), and modified-2 POX (MPX2), were applied to the OS part. The mutation operators for the MS part and the OS part were random replace and swap mutation, respectively.

Zhang et al. [25] used two-point crossover for the MS part and POX for the OS part. In order to maintain diversity in the population, the authors applied greedy mutation based on an SPT machine to the MS part. Moreover, they introduced a local search-based mutation operator called the adaptive neighbor search method for the mutation of the genes in the OS part.

Table 1 shows that the minimization of the makespan is the most popular objective function in the FJSP literature. This paper also considers minimization of the makespan of the FJSP. In addition, all of the previous research papers utilized a multichromosome that contains the MS and OS parts, which is also the case in this paper.

Most of the relevant research papers applied different crossover operators to the OS and MS parts of the FJSP multichromosome. Various previous research papers applied a GA with a single crossover operator to the FJSP [26–28]. However, they typically focused on one of the OS and MS parts, and details of the other part are missing. In such papers, the issue of genetic operators for the multichromosome was not discussed appropriately. In contrast, this paper applies a single crossover operator to the OS and MS parts in an evident manner.



Certain research papers developed GA hybrid solution methods and other algorithms for solving the FJSP, which are not considered in Table 1. Similar to the GAs in Table 1, almost all of the hybrid methods also adopted different genetic operators to handle the MS and OS parts.

The most important GA feature proposed in this paper is that crossover operations for both the MS and OS parts are performed by applying an identical operator: the candidate order-based genetic operator (COGO) provided by COGA. This can help to maintain a combination of a good MS gene and a good OS gene in a single operation during the search procedure. Moreover, we can see that the mutation for the OS part is also performed using the COGO. Note that the original version of the COGO is an integrated genetic operator, which can be used for the purposes of both crossover and mutation [14]. Consequently, the u-COGA for the FJSP can be developed by using only two kinds of genetic operators: the COGO and a mutation operator for the MS part. In other words, the structure of the u-COGO is simpler than that of previous GAs for the FJSP.

#### *2.2. Candidate Order-Based Genetic Algorithm*

The COGA is a type of GA that is based on the candidate order approach (COA). The main idea of the COA is to generate a feasible sequence of given items by appending one item at a time. Moreover, which item to append is chosen from a set of candidates; an item is a candidate if it is not appended to the sequence under construction and no constraints are violated after its appending [34]. COA can be used to develop metaheuristic algorithms for solving the sequencing problem.

The first example of a COA application is the COGA for solving the classical JSP (Kim, 2014). As a result of its flexibility, the COGA is quite effective for solving sequencing problems with additional constraints. Kim [14] applied the COGA to solve the job sequencing problem and the traveling salesman problem with precedence constraints. Kim and Kim [35] developed COGA for solving a variant of shortest path problem. Another important benefit of COA is that it can be integrated with various metaheuristic algorithms. For instance, Kim [34] developed a candidate order-based tabu search (COTS) for solving the job sequencing problem with precedence constraints. Kim and Kim [36] introduced the concept of the latest order constraint and applied COTS to solve the latest order constrained sequencing problem.

Typically, metaheuristic algorithms for solving sequencing problems have complex structures. However, COA-based metaheuristic algorithms are easier to implement than conventional metaheuristic algorithms. In particular, the COGA provides a distinguishing genetic operator called COGO, which has two important features. First, the COGO enables the application of simple point-based crossover strategies to sequencing problems, including one-point, two-point, and uniform strategies. Thus, both crossover operations for the MS and OS parts in a GA chromosome for the FJSP can be performed using the COGO. Second, the COGO also can be used to mutate a GA chromosome for the sequencing problem. The COGO generates a feasible sequence of given items in a constructive manner, which means a sequence is obtained by iteratively appending one item at a time. The item to be appended is chosen from a set of candidates. The objective of crossover, which is to create offspring similar to the parents, can be achieved by choosing an item with an earlier reference position. Note that the reference position of an item is its order in the parent solution. In contrast, the COGO performs the mutation by choosing an item with a later reference position, since the objective of the mutation is to create offspring dissimilar to the parents [14].

The main contributions of this paper are as follows. Firstly, we can implement the GA for solving the FJSP with only two genetic operators: the COGO and an additional operator for the mutation of the MS part. In other words, the u-COGA is easier to implement than previous GAs for the FJSP. Secondly, various kinds of crossover strategies, such as onepoint, two-point, and uniform strategies, can be applied to the u-COGA. Thirdly, numerical

experiments reveal that uniform crossover is the best crossover strategy for the u-COGA for the FJSP.

#### **3. Unified Candidate Order-Based Genetic Algorithm**

*3.1. Existing GAs for the FJSP*

Table 2 provides an example of the FJSP with three jobs (*J*1, *J*2, and *J*3) and three machines (*M*1, *M*2, and *M*3). Each job consists of two operations that should be processed in sequence, and *j*th operation of *Ji* (*i* = 1, 2, 3) is denoted by *oij*. An operation is processed by one of its alternative machines. For example, *o*<sup>11</sup> has three alternative machines (*M*1, *M*2, and *M*3). Different alternative machines can have different processing times for a single operation. In Table 2, *M*<sup>2</sup> is the SPT machine for *o*<sup>11</sup> such that its processing time is shorter than other machines. In contrast, *o*<sup>12</sup> has only two alternative machines (*M*<sup>2</sup> and *M*3), since *M*<sup>1</sup> is not available for this operation.


**Table 2.** An example of the FJSP.

In order to create a schedule for the FJSP, two kinds of decisions are needed. First, each operation should be assigned to one of its alternative machines. Second, operations should be sequenced. Typically, GAs for the FJSP utilize a multichromosome that consists of the OS part and the MS part for representing the operation sequence and machine assignment, respectively. An example of a multichromosome for the FJSP is shown in the upper part of Figure 1. Let *MS*- *oij* denote the machine to be used to process *oij*. It is clear that *MS*- *oij* must be an element of a set of alternative machines for *oij*.

**Figure 1.** Decoding the multichromosome representation for the FJSP.

In this paper, a sequence of operations indicates the assignment orders. Let *OS*- *oij* denote the order of *oij* in a sequence of given operations. Note that the encoding scheme for the OS part is the position listing representation, in which each gene specifies an order of the associated item. Thus, the OS part in Figure 1 can be converted into a sequence of given operations, *os*11, *os*12, *os*21, *os*31, *os*32, *os*22. An operation *oab* is scheduled prior to *ocd* if *OS*(*oab*) < *OS*(*ocd*). For instance, *o*<sup>11</sup> is the first operation to be scheduled in the Gantt chart in the lower part of Figure 1, since *OS*(*o*11) = 1. Moreover, we can start *o*<sup>11</sup> at time *t* = 0.

The order of given operations, *OS*- *oij* s, must satisfy the precedence relationships among operations within a single job. In other words, *OS*- *oij* must be smaller than *OS*(*oik*) if *j* < *k*. In Figure 1, *o*<sup>12</sup> is the second operation (*OS*(*o*12) = 2) of a job and its start time is 4, because we can start *o*<sup>12</sup> after its predecessor, *o*11, is completed. Valid *MS*- *oij* and *OS*- *oij* values yield a feasible schedule for the FJSP, which can be represented in the form of a Gantt chart, as shown in Figure 1.

The crossover operation is the "backbone" of the GA, and crossover operators are designed to generate two offspring solutions by recombining the genes of two parent solutions. The offspring solutions inherit the features of their parents. In other words, crossover operators should generate offspring solutions similar to their parents. However, this is difficult to achieve when solving the FJSP by applying a GA based on a multichromosome.

An example of a crossover operation of the existing GA for the FJSP is illustrated in Figure 2. Consider two parent solutions, *P*<sup>1</sup> and *P*2, in panel (a) in Figure 2. Typically, existing GAs for the FJSP apply different kinds of crossover operators to the OS and MS parts of parent solutions. In panel (b) in Figure 2, POX and classical one-point crossover are applied to the OS and MS part, respectively. POX is a well-known crossover operator for sequencing problems. In panel (b), the OS parts of the solutions are represented in the form of a sequence of operations, in order to apply POX. Job set 1 indicates that an operation *oij* has an identical position *OS*- *oij* in both *Px* and *Cx* (*x* = 1, 2) if it belongs to *J*3, where *Cx* denotes an offspring solution. In contrast, job set 2 means *C*<sup>1</sup> (*C*2) inherits precedence relationships among operations of *J*<sup>1</sup> and *J*<sup>2</sup> from *P*<sup>2</sup> (*P*1).




**Figure 2.** An example of crossover for an existing GA for solving the FJSP.

The one-point crossover operator divides a chromosome into two subparts, i.e., the earlier part and later part, and generates *C*<sup>1</sup> (*C*2) by combining the earlier part of *P*<sup>1</sup> (*P*2) and the later part of *P*<sup>2</sup> (*P*1). Panel (c) in Figure 2 shows a multichromosome representation of *C*<sup>1</sup> and *C*2, obtained by applying POX and one-point crossover. Similarly, most existing GAs

for solving the FJSP use different crossover operators to recombine the multichromosome, which consists of the OS and MS parts.

In Figure 2, we can see that the OS part and MS part evolve separately. That is, the crossover procedure of the OS part is not correlated with that of the MS part. In this context, the traditional approach based on two different crossover operators is called the separate evolution strategy in this paper. This strategy has a critical limitation, i.e., a combination of the OS gene and MS gene for an identical operation can be easily broken. For instance, *o*<sup>31</sup> has two associated genes: *OS*(*o*31) and *MS*(*o*31). For *P*<sup>1</sup> in panel (a) in Figure 2, *OS*(*o*31) = 1 and *MS*(*o*31) = *M*2. If these two genes indicate the optimal assignment order and machine for *o*31, the combination of *OS*(*o*31) = 1 and *MS*(*o*31) = *M*<sup>2</sup> should be preserved during the crossover procedure. However, this combination is not found in offspring solutions *C*<sup>1</sup> and *C*<sup>2</sup> in panel (c) in Figure 2. In more detail, *OS*(*o*31) of *P*<sup>1</sup> is inherited to *C*1, whereas *MS*(*o*31) of *P*<sup>1</sup> is inherited to *C*2. The u-COGA proposed in this paper is designed to overcome the limitation of the traditional separate evolution strategy for the FJSP.

#### *3.2. Crossover of Unified COGA*

This section outlines crossover procedure of the u-COGA. The COGO enables the application of simple crossover strategies to sequencing problems, such as one-point, twopoint, and uniform crossovers. Thus, a solution of the FJSP is converted into a sequence of operations in order to apply the COGO, as shown in Figure 3. Note that an operation in the FJSP has an additional attribute, *MS*- *oij* , as shown in panel (b) in Figure 3.



**Figure 3.** Reference set assignment for parent solutions in the form of a sequence.

The COGO generates a solution in a constructive manner, which means that one operation is appended to a new solution at a time. An operation to be appended to a specific position is chosen with consideration of the reference solution, one of the parent solutions. Thus, the COGO determines the reference solution for each position in a sequence before generating a new solution. An important feature of the COGO is that simple crossover strategies can be used to determine reference solutions. In panel (b) in Figure 3, each sequence is divided into two parts, i.e., the earlier part with the first two operations and later part with the other four operations, by adopting the conventional one-point crossover

strategy. We can see that *P*<sup>1</sup> (*P*2) and *P*<sup>2</sup> (*P*2) are used as reference solutions for the earlier part and later part of *C*<sup>1</sup> (*C*2), respectively.

An operation to be used as the *p*th position of a new solution is determined after the *p* − 1th operation is appended. An operation to be the *p*th operation is chosen from a set of candidates; an operation is a candidate if and only if it has not been appended to the new solution under construction and no constraints are violated by appending the operation to the *p*th position. The main idea of the COGO is that the objectives of crossover and mutation can be achieved by choosing the appropriate operation from the set of candidates.

Figure 4 illustrates a procedure for generating the earlier part of *C*1, an offspring of *P*<sup>1</sup> and *P*<sup>2</sup> in Figure 3. Let *RP*- *oij*, *X* denote the reference position of an operation *oij* in a solution *X*. In other words, *oij* is the *RP*- *oij*, *X* th operation in *X*. The main idea of the COGO is that the objectives of crossover and mutation can be achieved by choosing candidates with the smallest and the largest reference position values, respectively. For the first position in *C*1, all operations are candidates, since no operation is appended to *C*1, yet. Among the candidates, an operation with the smallest reference position regarding *P*1, *o*<sup>31</sup> is chosen as the first operation. Similarly, *o*<sup>11</sup> is appended to the second position in *C*1. Note that *MS*- *oij* s of the offspring are inherited from the associated reference solution. That is, *MS*(*o*31) and *MS*(*o*11) of *C*<sup>1</sup> are the same as those of *P*1.

**Figure 4.** Generating the earlier part of the offspring solution.

Figure 5 shows the crossover procedure for the later part of *C*1, where *P*<sup>2</sup> is used as the reference solution. We can see that the positions of *o*12, *o*21, *o*22, and *o*<sup>32</sup> in the later part of *C*<sup>1</sup> are different from those of *P*2. Nevertheless, the precedence relationships among them in *P*<sup>2</sup> are also maintained in *C*1. Moreover, *MS*- *oij* s of the operations in the later part are inherited from *P*2. Consequently, *C*<sup>1</sup> inherits from both *P*<sup>1</sup> and *P*2. Since the *MS*- *oij* value of the reference solution is maintained in the offspring solution, the crossover procedure described in Figures 3–5 helps to preserve a combination of the OS and MS genes.

In contrast, the mutation operation for the *p*th position of a new solution is performed by choosing a candidate with the largest reference position in the COGO procedure. That is, crossover and mutation are integrated in the COGO. Assume that *oij* is the candidate with the largest reference position for the *p*th position of a new solution chosen by the COGO in mutation mode. Then, *OS*- *oij* of the new solution is likely to be dissimilar to that of the reference solution. However, *MS*- *oij* is not affected by the mutation procedure provided by the COGO. Hence, we need an additional mutation operator for *MS*- *oij* , and


this paper applies random alternative mutation to the MS part of the multichromosome for the FJSP.

**Figure 5.** Generating the later part of the offspring solution.

#### *3.3. Initialization, Fitness Function, Selection, and Elitism*

Population initialization is the first step in the GA. Population can be defined as a set of chromosomes. There are several ways to initialize the population in the GA; however, the solutions in the initial population are randomly generated in this paper.

The fitness function of the GA is an objective function that needs to be maximized. In order to minimize the makespan of the FJSP, we use the reciprocal of the makespan as the fitness function of the u-COGA.

Selection is the procedure of creating a mating pool by copying the solutions in the current generation. The goal of the selection procedure is to choose the solutions with higher fitness values with larger probabilities. In this paper, we use conventional roulette wheel selection, in which a solution in the current population is chosen with a probability proportional to its fitness value.

Furthermore, the elitism strategy is applied to the u-COGA proposed in this paper. This strategy ensures that the best solutions in the current generation are also included in the mating pool. In other words, the elitism strategy helps to prevent desirable features of solutions being lost during the search procedure. Consequently, some solutions in the mating pool are determined by the elitism strategy, while the others are chosen by the roulette wheel selection procedure.

#### **4. Numerical Experiments**

In this paper, the u-COGA was written in Java language and was tested in a Windows 10 (64-bit) environment with an AMD Ryzen 7 2700X eight-Core processor (3.7 GHz) and 32 GB memory. The performance of the u-COGA was evaluated using the benchmark problems for the FJSP provided by the well-known Hurink library [37]. Depending on the average number of alternative machines for each operation, benchmark problems in the Hurink library are divided into three groups: eData, rData, and yData, as listed in Table 3.


**Table 3.** Groups of benchmark problems in the Hurink library.

Benchmark problems from each data group are summarized in Table 4, where *n* is the number of operations, *m* is the number of machines, and *Mij* is a set of alternative machines for an operation *oij*. Moreover, *Avg Mij* and *Max Mij* indicate the average and the maximum size of *Mij*, respectively. Each data group contains three instances, and the instances of the yData group have the most complex structures, in that they generally have *Avg Mij* and *Max Mij* larger than the other groups. On the contrary, eData contains the simplest instances.

**Table 4.** Benchmark problem instances.


As shown in Table 5, we applied three types of u-COGA, SS, TT, and UU to the benchmark problems in Table 4. The symbols 'S', 'T', and 'U' indicate one-point, two-point, and uniform crossover, respectively. In addition, the first and second symbols in 'type' indicate the crossover strategies for the OS part and MS part, respectively. For example, the first type of u-COGA, i.e., SS, applies the one-point crossover strategy to both the OS and MS parts.

**Table 5.** Types of crossover operators applied to the existing GA and u-COGA.


For comparison, the existing GA was also applied to the benchmark problems. The existing GA uses POX to recombine the OS parts of the chromosomes, while crossover for the MS part is performed in a one-point, two-point, and uniform manner. Thus, we have three types of existing GAs: PS, PT, and PU, as shown in Table 5.

The experiment results obtained by applying the u-COGA and the existing GA are summarized in Tables 6–8, where 20 experimental repetitions were performed for each case. *Cmax* and SD(*Cmax*) are the average makespan and standard deviation of the makespan, respectively. The six types of GAs applied to a single benchmark problem instance are ranked in ascending order of *Cmax*. For instance, TT and PS are the best and worst algorithms for the mt06 instance of eData in Table 6, respectively. All experimental results were obtained with a population size = 50, a crossover rate = 0.8, and a mutation rate = 0.01. In

addition, the elitism policy was used to preserve the five best solutions in the previous generation. A single experiment was terminated when its iteration number reached the maximum iteration limit = 300.


**Table 6.** Experiment results for eData.

**Table 7.** Experiment results for rData.



**Table 8.** Experiment results for vData.

Tables 6–8 provide the following observations. First, u-COGAs generally produce shorter makespan values than existing GAs, which indicates that the u-COGA is a promising approach for solving the FJSP. Second, for seven out of nine instances in Tables 6–8, two instances (mt06 and mt10) of eData, three instances (mt06, mt10, and mt20) of rData, and two instances (mt06 and mt10) of vData, the shortest makespan values were found by the UU-type u-COGA. In other words, uniform crossover was the best crossover strategy for minimizing the makespan of the FJSP. The performances of SS- and TT-type u-COGAs were not as good as that of the UU-type. Third, excepting one benchmark problem instance, i.e., mt06 of vData in Table 8, the worst makespan value for each instance was obtained by one of the existing GAs, which reveals the limitation of the separate evolution strategy for solving the FJSP. In particular, PS- and PU-type GAs produced the worst makespan values for four instances and three instances, respectively. Finally, the u-COGA and existing GA did not demonstrate significant differences in SD(*Cmax*) values.

Table 9 shows the average rank of the algorithms listed in Table 5. Again, we can see that the UU-type u-COGA has the smallest average rank, which means that it generally produces a shorter makespan than the other algorithms. In contrast, the average ranks of the individual existing GAs are larger than those of the u-COGAs. Consequently, the average rank of all u-COGAs (2.33) is noticeably smaller than that of all existing GAs (4.63), and we can conclude that the unified evolution strategy of the u-COGA helps to find better solutions for the FJSP.

**Table 9.** Average ranks of the algorithms.


Table 10 compares the best makespan of the existing GAs and u-COGAs. For instance, the best makespan for the existing GAs for mt20 of vData is 1068.85, while the best makespan for a u-COGA for the same instance is 1053.20, as shown in Table 8. Moreover, the reduction ratio in Table 10 was calculated as follows: (the best makespan of the existing Gas—the best makespan of u-COGAs)/(the best makespan of the existing GAs). This index was used to measure how much improvement was achieved by the u-COGA.



In Table 10, we can see that the u-COGA always produces better solutions for the FJSP than the existing GA. The reduction ratio values range from 1.36% to 5.29% (average reduction ratio = 2.79%). That is, the existing GA sometimes demonstrates a similar performance to the u-COGA; however, the u-COGA generally produces superior results for the FJSP.

#### **5. Conclusions and Further Remarks**

In order to apply the GA to solve a combinatorial optimization problem, a solution of the problem must be represented in the form of a chromosome comprising a number of genes. A solution for the classical JSP can be encoded using only OS-type genes. In contrast, a solution for the FJSP is typically represented in the form of a multichromosome, which consists of the OS and MS parts. Since the structures of the OS and MS parts are quite different, existing GAs generally apply different kinds of genetic operators to them. Consequently, the evolution procedures of the two parts are independent from each other, and a combination of OS and MS genes for an identical operation is easily broken during the search procedure. Moreover, the application of different kinds of crossover operators can make the GA complex and hard to implement and maintain. In order to overcome these problems, this paper proposes the u-COGA, which can be used to solve the FJSP.

The u-COGA utilizes an integrated genetic operator called the COGO to handle the multichromosome for the FJSP. The COGO enables the application of conventional crossover strategies, such as one-point, two-point, and uniform crossovers to the OS part. The u-COGA is designed to apply an identical crossover strategy to both the OS and MS parts. Moreover, the COGO generates a new sequence of given operations by choosing one operation at a time, which is suitable for the purposes of crossover and mutation. In order to perform the crossover operation, the COGO chooses a candidate with the smallest reference position. On the contrary, the mutation operation is performed by choosing a candidate with the largest reference position. In other words, the COGO of the u-COGA is used to perform three operations: crossover for the OS part, crossover for the MS part, and mutation for the OS part. Mutation for the MS part should be performed in a different manner, and random alternative mutation is used in this paper. Consequently, crossover and mutation are performed using two genetic operators: the COGO and mutation operator for the MS part in the u-COGA. Note that four different kinds of operators, i.e., crossover for the OS part, crossover for the MS part, mutation for the OS part, and mutation for the MS part, are required to implement the existing GAs for the FJSP. In other words, the COGO enables one to solve the FJSP using a GA with a simpler structure.

In order to evaluate performance of the u-COGA, numerical experiments were performed using an FJSP benchmark problem library. The experiment results reveal that the performance of the u-COGA is generally superior to that of the existing GAs in terms of the makespan. In more detail, the u-COGA improves the makespan values for the FJSP by about 2.79%. Notably, the best solution for a benchmark problem instance was generally found using the u-COGA based on the uniform crossover strategy. In contrast, the worst

solution was generally found by one of the existing GAs. This means that the u-COGA can explore the search space of the FJSP more effectively, and we can conclude that the GA for the FJSP should contain procedures for maintaining good combinations of good OS and MS genes; thus, the u-COGA is a promising approach for solving complex combinatorial optimization problems.

The authors propose three directions for future research on the u-COGA. First, this paper only considered one objective function, which is the minimization of the makespan. Thus, u-COGA performance regarding other objective functions should be investigated. Second, the primary advantage of the original version of the COGA is that it can flexibly handle additional constraints, such as precedence or position-related constraints. In this regard, in future research, we plan to apply the u-COGA to the FJSP with additional constraints. Third, it is expected that the u-COGA can help to solve other combinatorial optimization problems such that solutions have to be represented in the form of a multichromosome. Therefore, the authors will attempt to apply the u-COGA to combinatorial optimization problems in various other fields.

**Author Contributions:** Conceptualization, H.-Y.N. and J.-W.K.; methodology, J.-S.P. and J.-W.K.; software, J.-S.P.; validation, J.-S.P. and J.-W.K.; formal analysis, J.-S.P. and J.-W.K.; investigation, J.-S.P. and J.-W.K.; resource, H.-Y.N., T.-J.C. and Y.-T.N.; data curation, J.-S.P.; writing—original draft preparation, J.-S.P. and J.-W.K.; writing—review and editing, J.-S.P., H.-Y.N. and J.-W.K.; supervision, J.-W.K.; project administration, H.-Y.N., Y.-T.N. and J.-W.K.; funding acquisition, H.-Y.N., Y.-T.N. and J.-W.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (NRF-2020R1F1A1049251) and the Korea Institute for Advancement of Technology (KIAT) grant funded by the Korea Government (MOTIE: Ministry of Trade Industry and Energy; Grant No. N0002429).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Calibration of GA Parameters for Layout Design Optimization Problems Using Design of Experiments**

**Vladimir Modrak 1,\*, Ranjitharamasamy Sudhakara Pandian <sup>2</sup> and Pavol Semanco <sup>3</sup>**


**Abstract:** In manufacturing-cell-formation research, a major concern is to make groups of machines into machine cells and parts into part families. Extensive work has been carried out in this area using various models and techniques. Regarding these ideas, in this paper, experiments with varying parameters of the popular metaheuristic algorithm known as the genetic algorithm have been carried out with a bi-criteria objective function: the minimization of intercell moves and cell load variation. The probability of crossover (A), probability of mutation (B), and balance weight factor (C) are considered parameters for this study. The data sets used in this paper are taken from benchmarked literature in this field. The results are promising regarding determining the optimal combination of the genetic parameters for the machine-cell-formation problems considered in this study.

**Keywords:** facility layout; optimization; metaheuristic algorithm; cell formation; design of experiments

### **1. Introduction**

In general, facility-layout optimization problems are nonlinear, nonconvex, and multimodal in their nature. Facility-layout problems (FLP) can be divided according to types of manufacturing systems into four basic categories, which are product layout, process layout, static layout, and cellular layout [1]. Taking this classification into account, the proposed study addresses the cellular manufacturing problem. In the past, the main objective of a facility-layout problem was to minimize the material handling cost of the manufacturing system [2]. Presently, the main goal of the FLP is to improve manufacturing efficiency. In line with this, a number of authors suggested different objective functions for facility-layout problems, e.g., to maximize the throughput rate and minimize the conveyance time per trip [3] or minimize cycle times in order to increase productivity [4], and so forth.

An important issue in the design of cellular manufacturing systems (CMS) is the manufacturing-cell-formation problem (MCFP), which is based on group technology principles. Several taxonomies of the MCFP are reviewed in the literature, e.g., by Selim et al. [5], Bidanda et al. [6] and Modrak and Pandian [7]. The core of MCFP procedures is that machines are grouped into machine cells and parts into part families. In practical applications, it is not easy to arrange all parts and machines into autonomic cells, and therefore some operations have to be performed on separate machines. The cost of duplicating machines is often high, and therefore, related managerial decisions are usually trade-offs between economic and technological criteria [8–10]. During previous decades, numerous heuristic and metaheuristic methods and their variations were developed, tested and compared for this problem. The metaheuristics include the genetic algorithm (GA), simulated annealing, ant colony optimization, tabu search, scatter search, particle swarm optimization, GRASP and hybridized metaheuristics. There are several survey papers,

**Citation:** Modrak, V.; Pandian, R.S.; Semanco, P. Calibration of GA Parameters for Layout Design Optimization Problems Using Design of Experiments. *Appl. Sci.* **2021**, *11*, 6940. https://doi.org/10.3390/ app11156940

Academic Editor: Paolo Renna

Received: 3 July 2021 Accepted: 22 July 2021 Published: 28 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

such as by Herroelen et al. [11], Kolisch and Hartmann [12] and others, that present and compare the methods from various perspectives.

In this work, the focus is on application of GAs, and especially on finding out how the solutions of the cell-formation problem are influenced by a set of probability parameters of genetic operators, namely crossover and mutation, including balanced weight factors. In view of this, the presented work attempts to employ the Taguchi approach to find an optimal combination of parameters that impact the efficiency of the genetic algorithm and to explore whether the optimal combination of the genetic operators for the given type of MCFP can be influenced by the magnitude of the noise factors, which is represented by matrix size in our case.

#### **2. A Brief Literature Review**

In this section, a short review of the selected research on manufacturing cell formation will be given. John Holland, inspired by population genetics, introduced the concept of the GA in the 1970s. In 1975, the book [13] is published by him and his colleagues. The first GA-based approach for the cell-formation research area was proposed by Venugopal and Narendran [14]. It was a proposed mathematical model with a solution procedure based on a GA that can be purposely implemented in a cellular manufacturing environment. Joines and Houck [15] in their work applied a non-stationary penalty to solve CFP with a genetic algorithm. Gupta et al. [16] used a GA approach to solve the layout design problem with a predetermined number of manufacturing cells. Alsultan and Fedjki [17] utilized a GA approach to solve the machine-cell-part clustering problem in order to minimize total intercell and intracell moves. Gravel et al. [18] developed a GA with a double-loop, able to solve large-scale capacitated cell-formation problems with multiple routings. Moon and Gen [19] proposed a genetic algorithm to solve an integer programming model with consideration of alternative process plans and machine-duplication consideration. An adaptive genetic approach to solve the manufacturing-cell-formation problem in order to enhance the performance of the genetic search process was proposed by Mak et al. [20]. Zhao and Wu [21] evaluated the solutions of a multi-objective GA that applies minimizing costs due to intercell and intracell part flows, minimizing the total within-cell load variation and minimizing exceptional elements. They also incorporated in their research the multiple routes of parts. Arzi et al. [22] proposed a genetic algorithm with grouping efficiency and capacity requirements as objectives for large-scale systems design. Zolfaghari and Liang [23] tested and evaluated a genetic algorithm against simulated annealing and tabu search using binary cell-formation problems. Yasuda et al. [24] in their research proposed a method to solve the multi-objective CFP, partially adopting Falkenauer's grouping genetic algorithm. It was also aimed at improving the efficiency of their algorithm with regards to initialization of the population, fitness valuation, and keeping the crossover operator from cloning. Other similar approaches were published with many innovative algorithms. For example, Farahani et al. [25] described ant colony optimization to solve machine–part cell-formation problems. Mohammad Mohammadi et al. [26] approached a layout problem in cellular manufacturing systems with alternative processing routings. Dmytryshyn et al. [27] suggested a novel modeling approach for solving the cell-formation problem. Kamalakannan et al. [28] developed a simulated annealing algorithm for solving CFP with ratio level data. Shashikumar et al. [29] determined the solution for CFP using a heuristics approach. Sharma et al. [30] had done research on an implementation model using AHP and ANP for CFP. Octavio et al. [31] developed metaheuristic algorithms to solve grouping problems. Mourtzis et al. [32] brought the idea of adaptive scheduling in cellular manufacturing systems. Firouzian et al. [33] developed an artificial immune system for part family clustering. Vitayasak et al. [34] proposed a GA-based layout design approach to solve robust machine layout design problems for systems subject to demand uncertainties and maintenance. Their innovative method includes an experimental design that was used to test the robust design approach with corrective, preventative, and combined maintenance regimes. Dolgui et al. [35] explored a reconfigurable manufacturing system

that exhibited some crucial design and control characteristics for complex value-adding systems in highly dynamic scenarios. Such systems are designed at the outset for rapid change in the structure of machines, in order to quickly adjust production capacity and functionality within a part family in response to frequent market changes or intrinsic system changes. An interesting and very useful approach to the facility-layout design optimization was proposed by Bucki and Suchanek [36]. They proposed an effective tool for performance analysis of manufacturing systems from logistics viewpoint by using a mathematical simulation model.

It is worth mentioning other works that directly relate to the manufacturing cellformation problem, such as the hybrid GA/branch and bound approach to solve the manufacturing-cell-formation problem using a graph partitioning formulation, which was proposed by Boulif and Atif [37]. Their effort has been made to take into account the natural constraints of real-life production systems. A typical CFP with the objective of minimizing the exceptional elements was explored by Mahdavi et al. [38]. Deljoo et al. [39] used a GA to solve the dynamic cell-formation problem. Based on their studies, they reconsidered the shortcomings related to machine relocation cost and machine purchasing cost and developed a model for dynamic CFP in order to resolve these two shortcomings. Arkat et al. [40] developed a multi-objective GA for CFP considering cellular layout and operations scheduling. Cell-formation problems related to scheduling problems with the objective to minimize makespan are available in references [41–44]. A new algorithm for CFP with alternative machines and multiple-operation-type machines was developed by Li [45]. Its purpose was to improve traditional group technology cell-formation methods by considering alternative machines and multiple-operation-type machines. Boulif [46] proposed a new graph-cut-based encoding representation in order to solve the CFP with the genetic algorithm. Obviously, there are other related works, since this domain attracts a large research interest.

#### **3. Structure of Genetic Algorithm**

In a GA, a high-quality candidate solution is represented by a collection of genes called chromosome. A chromosome's potential is given by its fitness function. A population consists of a set of selected chromosomes, and the population is subjected to generations (or iterations). Finally, crossover and mutation operators are performed with defined probabilities to improve the solutions. A GA has several advantages over the traditional optimization methods. It may quickly arrive at a good solution set. As worse cases are eliminated, they will never affect the generated solution. GAs are one of the best methods to solve a problem about which little is known. This is because a genetic algorithm works by its own rules. This is a very useful strategy for a GA to solve highly complex problems in nature. In a GA, it is necessary for the problem solver to choose the appropriate coding method. If the solutions are coded in different combinations, then the GA will start its searching operation using its operators known as selection, crossover, and mutation, respectively. As is known, a suitable type of crossover technique for a particular problem can improve the GA's performance [47]. All these methods are probabilistic in nature. The proper stopping criteria will be given as input for the GA to stop its searching process. This is done purely based on the experience of the problem-solver. Based on the stopping criteria, the GA will stop running and give the solution that it finds at that point of time. The solution of the problem has to be represented in the GA as a genome (or chromosome). The genetic algorithm then creates a population of solutions and applies genetic operators such as mutation, crossover, and selection to evolve the solutions in order to find the best one(s). The crucial aspects of using GAs are: the definition of the objective function, the definition and implementation of the genetic representation, the definition and implementation of the genetic operators [48].

GAs are frequently adopted to find out how to make the machine clusters form cells in accordance with the rules of production flow analysis. In cell formation, a GA works well and finds a good solution, as given in the abovementioned literature studies. It gives high-quality solutions even if the problem has a high level of complexity. There are few other traditional metaheuristics approaches that performing good searches in such a situation [49]. The following representation is used in a typical cell-formation problem solved using GA. This representation is popularly known as real coding (see Figure 1).



**Figure 1.** Representation of a chromosome using real coding.

Depending on the number of cells to be accommodated in the layout, the number of genes on a chromosome increases. In the above example there are three cells and hence there are only 1, 2, and 3 values present in the chromosome of 8 genes (8 machines).

#### **4. Mathematical Model of the Cell-Formation Problem**

According to the review of the literature, the minimization of intercell flows and the total cell load variation can be considered the essential objectives in the manufacturing cell-formation research. Thus, the bi-objective fitness function used to evaluate the solution incorporates intercell flows and cell load variation. The mathematical model is given as:

Intercell flow fitness function:

$$\text{Minimize}\\f\_1 = \sum\_{i} \sum\_{j} \sum\_{k} a\_{ij} |X\_{ik} - Y\_{ik}| \tag{1}$$

Cell load variation fitness function: Minimize

$$f\_2 = \sum\_{i} \sum\_{k} X\_{ik} \sum\_{j} \left(\omega\_{ij} - m\_{ik}\right)^2 \tag{2}$$

Variable:

$$m\_{\rm ik} = \sum\_{i} X\_{\rm ik} \times Y\_{\rm jk} \times \omega\_{\rm ij} / \sum\_{i} X\_{\rm ik} \tag{3}$$

Decision variables:

$$
\omega\_{\rm ij} = t\_{\rm ij}; \text{ if } \text{part } j \text{ is processed on machine } i \tag{4}
$$

0; otherwise

$$X\_{ik} = 1; \text{ if machine } i \text{ is in cell } k \tag{5}$$

*Yjk* = 1; if part *j* is in cell *k* (6)

0; otherwise

0; otherwise

$$a\_{ij} = 1; \text{ if } \text{part} \neq \text{needs to be processed on machine } i \tag{7}$$

0; otherwise Constraints: *<sup>c</sup>*

$$\sum\_{k}^{\mathcal{L}} X\_{ik} = 1 \; \forall \in \{1, 2, \dots, m\} \tag{8}$$

$$\sum\_{k}^{c} \mathbf{Y}\_{jk} = 1 \; \forall \; \in \; \{1, 2, \dots, p\} \tag{9}$$

$$\sum X\_{ik} \ge 1 \; \forall \; k \in \{1, 2, \dots, c\} \tag{10}$$

$$\sum Y\_{\vec{j}k} \ge 1 \,\,\forall \, k \in \{1, 2, \dots, c\} \tag{11}$$

Bi-objective fitness function:

$$Z(t) = \mathfrak{a} \cdot f\_1 + (1 - \mathfrak{a}) \cdot f\_2 \tag{12}$$

The given model of cell-formation problem works with processing time. According to this model, we search to obtain the number of cells and the number of machines and parts within each cell. Equation (1) shows the calculation of the intercell flows, and Equation (2) shows the cell load variation. Equation (3) shows the average intercell processing time for the *j*-th part and *k*-th cell. Equation (4) shows the decision variables that assign the time *tij* of the *i*-th machine required to process the *j*-th part. Equations (5)–(7) are decision variables that state that *xi*k, *yj*<sup>k</sup> and *ajj* are 0–1 binary numbers. Equations (8) and (9) ensure that each machine and part is attached to only one cell. Equations (10) and (11) ensures that, in each cell, there must be allocated at least one machine and one part, respectively. Equation (12) is the bi-objective fitness function for a non-binary cell-formation problem balanced by weight factor *α*.

A fitness function value is computed for each chromosome in the population, and the objective is to find a chromosome with the maximum fitness function value. Due to objective of minimizing both the total cell load variation and the exceptional elements, it is necessary to map it inversely and then maximize the result. Goldberg [50] suggested a mapping function given as:

$$F(t) = Z\_{\text{max}} - Z(t) \tag{13}$$

The symbol *F*(*t*) stands for the fitness function of the *t*-th chromosome, and *Zmax* is the *max*[*Z*(*t*)] of all chromosomes (*t*). The advantage is that the worst chromosomes obtain a zero-fitness function value, so they are not going to be reproduced into the next generation.

The following procedure given by Zolfaghari and Liang [51] is used to assign parts into the machine cells:

$$P\_{kj} = \left(\frac{f\_{kj}}{f\_k}\right) \cdot \left(\frac{f\_{kj}}{f\_j}\right) \cdot \left(\frac{T\_{kj}}{T\_j}\right) \tag{14}$$

where, *Pkj* is the membership index of the *j*-th part belonging to the *k*-th cell; *fkj* is the number of machines in the *k*-th cell required by the *j*-th part; *fk* is the total number of machines in the *k*-th cell; *fj* is total number of machines required by the *j*-th part; *Tkj* is the processing time of the *j*-th part in the *k*-th cell; and *Tj* is the total processing time required by the *j*-th part.

#### **5. Case Study on Layout Design Optimization**

The procedure we decided to apply to analyzing the influence of genetic parameters on the final solution quality for reorganization of machines and parts into cells is an experimental design method developed by Genichi Taguchi. Here we established a Pdiagram (see Figure 2) that identifies the inputs and outputs of the system together with control and noise factors.

The Taguchi experimental method consists of performing selected experiments to study the influence of several operating factors on output-parameter values. Taguchi separates the factors into two domains: control and noise factors. The difference between these two factor groups is that we cannot control them directly. The elimination of the noise factors is often impractical and impossible so the Taguchi method seeks to minimize the effect of noise, using optimal levels of important control factors based on the concept of robustness [52]. These three fundamental control factors are the probability of crossover

(A), probability of mutation (B), and balance weight factor (C). The levels of each set of control factors are shown in Table 1.

**Figure 2.** Parameter diagram for a genetic algorithm.

**Table 1.** Input signal factors and levels.


Based on the principles of the genetic algorithm, these control-factor values are decided taking a reference from the optimization of control parameters for a genetic algorithms using an image registration problem [53]. In the evaluation of a final feasible solution, there also exist the so-called noise factors that have an impact on the results. These factors cause deviation in the search space size and the noise factors, such as the size of the machine-part matrix (D). The offered factors are composed of individual levels, wherein each level is completely independent. This particular case is under consideration with three levels of control and two levels of noise factors.

It is proposed to apply Taguchi quality concept for the assessment of the final solution, and, in this context, the L9 orthogonal array has been chosen. The results of the experiments are subsequently transformed into a signal-to-noise ratio. For reducing the variability of solutions around a target, the smaller-is-better S/N ratio (SNR) calculation is applied [54]:

$$SNR = -10\log[\frac{i}{n}\sum\_{i=1}^{n}Y\_i^2] \tag{15}$$

A signal-to-noise ratio is a measure of robustness that can be used to identify the control factor settings that minimize the effect of noise on the response. In the parameter design, there are two factors: the signal factor that can be controlled and the noise factor that is too expensive to control. In this research work, the signal factors are A, B, and C, wherein the noise factor is the exceptional element. *Yi*<sup>2</sup> is obtained from the mean sqare deviation (MSD) of signal factors and of exceptional elements. Here the '*n*' represents the number of experiments that are performed using Taguchi's experiment design.

Higher values of SNR will control the noise factors that lead to the minimization of the objective function. Since, in this research work, the objective is to minimize the combination of two functions, it is recommened from Taguchi's SNR principle that the 'smaller is better' is suitable for this approach.

The stopping criterion used to test algorithms was set at a number of generations [55], fixed to min 120. In order to conduct the experiment, we used a set of four instances and implemented the algorithms in PHP script. As a common performance measure, the number of exceptional elements (EE) has been used. The total set of experiments that were performed was obtained by combining the L9 array of the control factors. In order to obtain more objective results, we decided on two groups of size problems; the first of them consisted of 4 size problems (24 × 16, 24 × 40, 30 × 16, 30 × 40) as shown in Table 2, and the second one consisted of 19 size problems (see Table 3).


**Table 2.** Combination of factors and the resulting values of the trials of experiments (Group #1).

**Table 3.** Combination of factors and resulting values of the trials of experiments (Group #2).


The total number of experiments for the first group of size problems is 36 and for the second group of size problems is 171. The objective in this factorial experiment is to search values of the signal-to-noise ratio for "smaller-is-better", which is computed for each set of experiments. After obtaining the results of Taguchi's experiment design, EEs are transformed into S/N ratios. The results of the experiment for each series of experiments are presented in Tables 2 and 3.

#### **6. Conclusions**

Based on Table 1 (input signal factors and levels), the crossover probability is considered in ascending order (lower to higher), which is mentioned in Tables 2 and 3. The change in value for A occurs in every three rows; for B, the change in value occurs in every row in ascending order, and for factor C (1,2,3 then 2,3,1 then 3,1,2) as per the L9 orthogonal array representation. It happens incidentally that the last value (9th row) is better compared to other values, as far as this research concerned. The optimal level of the factors is the level with the highest SNR. Tables 2 and 3 show that the optimal level has been obtained in the 9th row (refer to Tables 2 and 3) where the value of factor A is 0.8, the value of factor B is 0.1, and the value of factor C is 0.6, mentioned as (+, +, 0) as per Table 1. The general view is that whenever the probability of crossover and/or mutation is increased, the number of searches in the solution space will also increase; thereby the chance of obtaining a better solution is increased. The Taguchi's design of experiments based on this research work provides evidence to the above statement.

Building on this analysis the following conclusions and suggestions for further research have been made.

It can be stated that the Taguchi design of experiments method can be used as an effective tool to determine the optimal combination of the genetic operators for the given type of MCF problems, because two different experiments brought out that the optimal level of the factors A, B, and C are identical.

Due to the limited number of experimental groups, we can only anticipate that an optimal combination of the genetic operators for the given type of MCF problem is not influenced by the magnitude of the noise factors. Therefore, in our further research this dependence/independence relation will be investigated.

This work provides further evidence that the efficiency of a genetic algorithm is dependent on the mentioned control factors and their parameters. As a direction for future studies, it could be interesting to extend the parameters (for example, the probability of reproduction and the number of generations) and develop effective genetic algorithm incorporating advanced features. For more realistic models, it would be useful to consider several practical assumptions, such as machine-availability constraints or changeover times.

**Author Contributions:** Writing—original draft preparation was done by R.S.P. Conceptualization and methodology were performed by all the authors. Writing—review was done by V.M. Experimentation was carried out by P.S. Results of the investigations were analyzed and evaluated by all the authors. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors would like to thank the Ministry of Education, Science, Research and Sport of the Slovak Republic for the financial support of this work by grant KEGA, project No. 025TUKE-4/2020.

**Data Availability Statement:** The detailed input numerical data used in this research work are available from the corresponding author upon request. Due to the actual fact that data were collected from various literature, it is used for research purposes only.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **An Algorithm for Arranging Operators to Balance Assembly Lines and Reduce Operator Training Time**

**Ming-Liang Li**

Department of Industrial Management, Asia Eastern University of Science and Technology, Banqiao Dist., New Taipei City 22061, Taiwan; fg017@mail.aeust.edu.tw; Tel.: +886-(2)-7738-0145 (ext. 5113)

**Abstract:** Industry 4.0 is transforming how costs, including labor costs, are managed in manufacturing and remanufacturing systems. Managers must balance assembly lines and reduce the training time of workstation operators to achieve sustainable operations. This study's originality lies in its use of an algorithm to balance an assembly line by matching operators to workstations so that the line's workstations achieve the same targeted output rates. First, the maximum output rate of the assembly line is found, and then the number of operators needed at each workstation is determined. Training time is reduced by matching operators' training and skills to workstations' skill requirements. The study obtains a robust, cluster algorithm based on the concept of group technology, then forms operator skill cells and determines operator families. Four numerical examples are presented to demonstrate the algorithm's implementation. The proposed algorithm can solve the problem of arranging operators to balance assembly lines. Managers can also solve the problem of worker absences by assigning more than one operator with the required skillset to each workstation and rearranging them as needed.

**Keywords:** assembly line balancing; group technology; cluster algorithm; bottleneck station; output rate

#### **1. Introduction**

Netessine and Taylor [1] demonstrated that expensive production technology usually leads to low product prices and may simultaneously lead to high-quality products. However, the product price is usually proportional to the product cost, part of which is labor cost. Assembly line balancing (ALB) is a crucial part of production/operations management [2]. On most mass-production assembly lines, workers repetitively perform a set of tasks predefined by using assembly line balancing techniques [3], which assign a number of work elements to various workstations to maximize the assembly line's balancing efficiency [4]. A balancing of tasks across the workstations results in the optimization of a given objective function without violating preceding constraints [5]. An unbalanced assembly line can be dangerous. For example, an operator without adequate training or who is focused on personal problems can cause an accident, decrease the efficiency of the assembly line, and increase production costs. Assembly lines maximize assembly operations, and most manufacturing plants have one or more lines, making assembly line balancing one of the oldest challenges in industry [6]. Moreover, a balanced assembly line can reduce a factory's work-in-process (WIP) inventory. Kim et al. [7] suggested a mathematically mixed integer programming model that accounts for the WIP balance and setup time in a semiconductor fabrication line. Maximum productivity and minimum load smoothness are benefits of solving the assembly line balancing problem (ALBP) [8,9].

Many types of production lines exist. Kuroda [10] presented a robust design for a cellular-line production system with unreliable facilities. A cellular-line production system comprises multiple flow shops, which are individual sets composed of functionally different facilities. Kuroda grouped facilities that perform similar operations.

**Citation:** Li, M.-L. An Algorithm for Arranging Operators to Balance Assembly Lines and Reduce Operator Training Time. *Appl. Sci.* **2021**, *11*, 8544. https://doi.org/10.3390/ app11188544

Academic Editors: Vladimir Modrak and Zuzana Soltysova

Received: 26 July 2021 Accepted: 12 September 2021 Published: 14 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

This study initially considers a traditional assembly line that has four workstations with outputs *λA*, *λB*, *λC*, and *λ<sup>D</sup>* (Figure 1). Thus, the production rate, or output rate, of this assembly line is equal to *λline* = *min*.{*λA*, *λB*, *λC*, *λD*}.

**Figure 1.** Traditional assembly line.

The cycle time (CT) refers to the longest period required for workstations to complete operations. Therefore, all station operation times are less than or equal to the CT. The aim of balancing the assembly line is to ensure that the output rate and operating time of each workstation are identical. Previous studies have examined assembly line balancing. According to Bartholdi and Eisenstein [11], the assembly line balances itself. They found that if workers are sequenced from slowest to fastest, a stable partition of work spontaneously emerges independent of the stations at which the workers begin. Moreover, the production rate converges to a value. Huang [12] determined that the CT and labor efficiency are negatively exponential to the number of assigned workers. Sawik [13] presented new mixed integer programming formulations for scheduling a flexible flow line with blocking. In this flexible flow line, each stage has one or more identical parallel machines. This study uses parallel operators instead of parallel machines. Rearranging operators may increase the production rate. Muth [14] examined whether a rearrangement in the order of servers affects the production rate. Muth speculated that the production rate remains invariant if the production line is reversed.

If the number of workers at a workstation is more than the number of machines, the machines are expected to limit the output rate and cause a bottleneck. However, Zavadlav et al. [15] showed that the workers, rather than the machines, limit the rate of output.

The manager of a factory always works under "cost" pressure; the company wants the manager to make continual improvements without increasing the cost. Managers typically try to increase the output by adding machines and workers, but managers must first consider whether assembly lines are balanced. When assembly lines are not balanced, productivity is low and workers complain. Workers might not only complain that they have to do more work at their station for the same pay, but also that their skills are not suited to their workstation. This work environment might lead workers to miss shifts on purpose or to resign.

Many ALBP methods have been examined in the study of manufacturing systems, including the heuristic algorithm [16–24], artificial bee colony algorithm [25–27], genetic algorithm [28–30], branch-and-bound algorithm [31–34], imperialist competitive algorithm [35], Pareto Greedy algorithm [36], memetic ant algorithm [37], whale optimization algorithm [38], recursive algorithm [39], fuzzy linear programming [40], swarm algorithm [41,42], simulated annealing algorithm [43], and data envelopment analysis [44]. The genetic algorithm [45,46] and artificial bee colony algorithm [47] have been used to solve the ALBP in remanufacturing systems. However, studies have seldom discussed how to arrange operators and group similarly skilled operators at the same workstation. This paper considers workers, not machines, in determining whether a workstation has the right equipment and machinery for operators to use. To achieve this goal, the study works out a method for rapidly balancing assembly lines by arranging operators among the workstations and uses a robust, cluster algorithm based on the concept of group technology (GT) to form operator skill cells and determine operator families.

#### **2. Problem Statement**

Traditionally, when arranging a limited number of operators on an assembly line, a manager begins with one operator per workstation, before adding an operator to the workstation with the lowest output. Step by step, all of the operators are arranged and the output rate can be found.

When a workstation has a lower output rate, the manager typically adds another operator. However, the lower output might be due to the "human" factor; an operator who is trained on many machines might be assigned to a workstation outside of their skill set. For example, if an operator knows how to work an auto-optical inspection machine but is assigned to a workstation with an in-circuit test machine, the station's output rate will drop because the operator does not have the required skills. Output can also fall when a worker is absent or resigns unexpectedly and the manager fills the vacancy with a random worker who has not been trained to operate that workstation.

Thus, the problem statement of this study is:

Ideal: The maximum output rate is found, all of the workstations have the same output rate, and the number of operators at each workstation can be found, with all operators assigned to a workstation that matches their skill set.

Reality: The arrangement algorithm is step by step, and the maximum output rate cannot be determined until all of the operators have been arranged. Operators can be assigned to unsuitable workstations.

Consequences: The maximum output rate is an unknown that will affect a boss's decisions, especially before setting up a new plant. An operator who works at an unsuitable workstation might complain, be absent, or quit unexpectedly, causing a lower output rate.

Proposal: An arrangement algorithm with a constant number of operators must consider the maximum output rate and balance the assembly line without any increase in cost. A sorting algorithm must consider the skills required at each workstation and the skill set of each operator.

#### **3. Methodology**

A manager adds a worker to a workstation that has a lower output than at other workstations. When the output of another workstation is lower, the manager adds another worker to that workstation. The use of this step-by-step method balances the assembly lines, especially when setting up a new plant.

#### *3.1. Example of Step-by-Step Method*

Suppose an assembly line has six workstations, as shown in Figure 2.

**Figure 2.** Assembly line with six workstations.

In addition, suppose that the terms *t*11, *t*12, *t*22, *t*13, *t*14, and *t*<sup>15</sup> are equal to 12, 15, 10, 16, 20, and 8 min, respectively. If one workday comprises 8 h, then the output per operator per day at each workstation, *γ*11, *γ*12, *γ*22, *γ*13, *γ*14, and *γ*15, are 40, 32, 48, 30, 24, and 60, respectively. Suppose that 500 workers are assigned to this assembly line. The workers can be assigned step by step, as follows (Table 1).

**Table 1.** Traditional ALB assignment.


The term *λline* is limited in workstation *a*14, which is a bottleneck workstation. Thus, this workstation should have one more operator to increase the output rate. The seventh worker is added to workstation *a*14, and the following is obtained (Table 2).

**Table 2.** Traditional ALB assignment.


The assignment of workers in a step-by-step manner uses considerable time to obtain the following final result (Table 3).

**Table 3.** Traditional ALB assignment.


Considerable time was spent and many steps performed to determine the optimal arrangement of operators at the assembly line's workstations. Although some studies have discussed and tried to solve the ALBP, the question remains: Is there another method for arranging the workers, instead of using the step-by-step method? This paper presents a novel approach for solving this problem.

Almost all workers require training before operating a workstation. Ideally, the operators on an assembly line would have all of the skills required for each workstation, but that would be very costly. Thus, the operators' skills must be determined so that they can be placed at a workstation that requires their skills. For example, consider an assembly line that has welding and drilling workstations. If operator κ possesses welding skills, assigning operator κ to the welding workstation rather than to the drilling workstation will reduce the training time for the welding workstation.

Suppose that there is a balanced assembly line with seven workstations (Figure 3), and workstations 1–7 require one, three, four, two, two, two, and one operator, respectively.

**Figure 3.** Assembly line with seven workstations.

Each workstation requires a skill, which is denoted as *s*1*, s*2, ..., *s*7. A matrix *SO* = *soij uv* describes the relationship between skills and operators. For example, *soij* = 1 indicates that operator *j* possesses skill *i*. For the assembly line in Figure 3 to function, the minimum matrix *SO* must be as follows:

ଵ ଽ ଼ ହ ସ ଷ ଶ ଵ ଵ ଵ ଵ ଶ ଵ ଷ ଵ ସ ଵ <sup>ହ</sup> ܱܵ ൌ ଵݏ ଶݏ ଷݏ ସݏ ହݏ ݏ ۏ ݏ ێ ێ ێ ێ ێ ۍ ͳ ͳͳͳ ͳͳͳͳ ͳ ͳ ͳ ͳ ͳ ͳ ے ͳ ۑ ۑ ۑ ۑ ۑ ې

Suppose it is known that:


On the basis of the first-come, first-served rule, operator 1 is assigned to *s*1; operators 2, 3, and 4 are assigned to *s*2; operators 5, 6, 7, and 8 are assigned to *s*3; operators 9 and 10 are assigned to *s*4; operators 11 and 12 are assigned to *s*5; operators 13 and 14 are assigned to *s*6; and operator 15 is assigned to *s*7. Thus, at least 11 operators must take the training course, namely operators 1, 3, 4, 5, 6, 8, 9, 11, 12, 14, and 15, as illustrated in the following matrix where an "X" denotes that the operator needs the training course:

$$SO = \s\_4 \begin{bmatrix} \underbrace{\mathbf{X}}\_1 & \mathbf{o}\_2 & \mathbf{o}\_3 & \mathbf{o}\_4 & \mathbf{o}\_5 & \mathbf{o}\_6 & \mathbf{o}\_9 & \mathbf{o}\_9 & \mathbf{o}\_1 & \mathbf{o}\_1 & \mathbf{o}\_1 & \mathbf{o}\_1 & \mathbf{o}\_1 & \mathbf{o}\_1 \\ \underbrace{\mathbf{s}\_1}\_3 & \mathbf{1} & \underbrace{\mathbf{1} & \mathbf{X} & \mathbf{X}}\_1 & & & & & \\ \mathbf{s}\_3 & & & & & & & & \\ \mathbf{s}\_5 & & & & & & & \\ & \mathbf{s}\_6 & & & & & \\ & & & & & & \\ & & & & & & 1 & 1 & 1 \end{bmatrix} \underbrace{\begin{bmatrix} 1 & \mathbf{1} & \mathbf{1} & \mathbf{1} & \mathbf{1} & \mathbf{1} \\ \mathbf{X} & \mathbf{X} & \mathbf{1} & \mathbf{X} \end{bmatrix}}\_{\mathbf{S}\_1} \underbrace{\begin{bmatrix} \mathbf{1} & \mathbf{1} & \mathbf{1} \\ \mathbf{1} & \mathbf{1} & \mathbf{1} \end{bmatrix}}\_{\mathbf{T}\_2} \underbrace{\begin{bmatrix} \mathbf{1} & \mathbf{1} & \mathbf{1} \\ \mathbf{1} & \mathbf{1} & \mathbf{1} \end{bmatrix}}\_{\mathbf{T}\_3} \underbrace{\begin{bmatrix} \mathbf{1} & \mathbf{X} \\ \mathbf{1} & \mathbf{X} \end{bmatrix}}\_{\mathbf{T}\_4} \underbrace{\begin{bmatrix} \mathbf{1} & \mathbf{X} \\ \mathbf{1} & \mathbf{X} \end{bmatrix}}\_{\mathbf{T}\_5}$$

A cluster algorithm is introduced in this paper to reduce the training time and increase the speed at which an assembly line functions. Workers' absences are also a problem for managers [48,49], but the proposed cluster algorithm can suggest which worker to rotate to a workstation with an absent operator.

Instead of using the traditional step-by-step method, this study aims to use the proposed algorithm to assign operators to each workstation. Before assigning the operators, some basic data must first be known about each workstation, such as the required skills, the standardization time, the work time per day, and the total number of operators. Then, the approximate number of operators needed at each workstation can be found by using Equations (1)–(7) in Section 3.1. Finally, the number of operators at each workstation can be adjusted by using the traditional step-by-step method. The training time can be reduced by assigning an operator with suitable skills to a suitable workstation, or by sorting the skill–operator matrix by using the cluster algorithm proposed in Section 3.2. The concept is shown in Figure 4.

**Figure 4.** The novel algorithm proposed in this study.

*3.2. Number of Operators at Each Workstation Decided*

Suppose an assembly line has *m* × *n* workstations (Figure 5).

**Figure 5.** An assembly line with *m* × *n* workstations.

The matrix *A* = *aij mn* can be used instead of the previous assembly line.

$$A = \left[ \begin{array}{cccc} a\_{11} & a\_{12} & \dots & a\_{1n} \\ & \vdots & \vdots & \ddots & \vdots \\ & a\_{m1} & a\_{m2} & \dots & a\_{mn} \end{array} \right]$$

The output rate of the assembly line is *λline* where

$$
\lambda\_{Iine} = \min \{ \lambda\_{i\bar{j}} \}, \ 1 \le i \le m, \ 1 \le j \le n
$$

The ideal assembly line means that

$$
\lambda\_{line} = \lambda\_{11} = \lambda\_{21} = \dots = \lambda\_{mn} \tag{1}
$$

because

$$
\lambda\_{i\dot{j}} = p\_{i\dot{j}} \gamma\_{i\dot{j}} \tag{2}
$$

The term *pij* can be easily determined, as follows:

$$p\_{ij} = \frac{\lambda\_{ij}}{\gamma\_{ij}} \tag{3}$$

Moreover, the total number of operators is *P*. , where

$$P = \sum\_{i=1}^{m} \sum\_{j=1}^{n} p\_{ij} \tag{4}$$

Thus, the following equations are obtained:

$$P = \sum\_{i=1}^{m} \sum\_{j=1}^{n} \frac{\lambda\_{ij}}{\gamma\_{ij}} = \lambda\_{lim} \cdot \sum\_{i=1}^{m} \sum\_{j=1}^{n} \frac{1}{\gamma\_{ij}} = \lambda\_{lim} \cdot \frac{1}{T} \sum\_{i=1}^{m} \sum\_{j=1}^{n} t\_{ij} \tag{5}$$

and,

$$\lambda\_{Iinc} = \frac{P}{\sum\_{i=1}^{m} \sum\_{j=1}^{n} \frac{1}{\gamma\_{ij}}} = \frac{P \cdot T}{\sum\_{i=1}^{m} \sum\_{j=1}^{n} t\_{ij}} \tag{6}$$

Replace *λij* from Equation (3) with *λline* from Equation (6) to obtain the following equation:

$$p\_{ij} = \frac{\lambda\_{ij}}{\gamma\_{ij}} = \frac{P}{\gamma\_{ij} \sum\_{i=1}^{m} \sum\_{j=1}^{n} \frac{1}{\gamma\_{ij}}} = \frac{P \cdot T}{\gamma\_{ij} \sum\_{i=1}^{m} \sum\_{j=1}^{n} t\_{ij}} \tag{7}$$

As term *pij* is not always an integer, it may need to be chopped off.

Almost all of the operators are now assigned. The traditional step-by-step method is adjusted (as demonstrated in the previous section) and used to assign the remaining workers.

#### *3.3. Cluster Algorithm*

The proposed cluster algorithm is based on the concept of GT, which has been applied to cellular manufacturing systems [50,51]. The concept involves grouping machines into machine cells and parts into part families [52–59]. Several algorithms, such as heuristic algorithms [60,61], genetic algorithms [62,63], and neural network algorithms [64], have been derived for solving the GT problem. The proposed algorithm begins with the closeness matrix and a comparison with other algorithms in example 3.

First, the relationship between *u* skills and *v* operators, which is represented by the matrix *SO* = *soij uv* is determined as follows:

$$so\_{ij} = \begin{cases} \begin{array}{c} 1, \text{ } if \stackrel{\circ}{\ } operator \ j \text{ has } i \text{ } skill \end{array} \\ \begin{array}{c} 0, \text{ } otherwise \end{array} \end{array} $$

The closeness matrix of skills is defined as *Bu* <sup>=</sup> {*brs*}*u*×*u*, where

$$b\_{rs} = \begin{cases} \begin{array}{c} \sum\limits\_{k=1}^{v} \text{so}\_{rk} \text{so}\_{sk} \quad \text{if } r \neq s\\ 0 & \text{, otherwise.} \end{array} \end{cases} \tag{8}$$

For example, if there are three workers that know how to operate the welding and milling machines, the closeness between the welding and milling skills is three.

#### 3.3.1. Skills Sorting

Whether two skills are grouped depends on their closeness. Closeness denotes how many operators have both of the skills. The more operators have the two skills, the closer the skills are. First, the skill with the highest closeness value is found and a new group is formed. Second, an unsorted skill follows a sorted skill, based on the highest value in its row in the closeness matrix. If there is no sorted skill to follow, the skill is placed in a new group. Finally, all skills are grouped.

Thus, the skills can be sorted in three steps:

Step 1. Let *f* ∈ *G*1, where *f* = { *x*|*bxk* = *max*.{*brs*}, 1 ≤ *k* ≤ *v*}. A tie is broken by selecting the smaller one.

Step 2. Ignore all the grouped rows and form a new closeness matrix *B u* = *bij* . Let

$$\begin{cases} \quad f \in G\_{k\prime} & \text{if } = \{ b\_{\vec{\imath}\vec{\jmath}} = \max \{ b\_{\vec{\imath}\vec{\jmath}} \} \land y \in G\_k\\ \quad h \in G\_{\text{new\\_group}\prime} & \text{otherwise} \end{cases}$$

Step 3. Repeat Step 2 until *u* skills are grouped.

#### 3.3.2. Operator Sorting

Operator sorting aims to place operators in the group that will use the highest number of their skills.

$$\begin{aligned} \text{Thus, let } v \in G\_{y, \epsilon} \ y &= \left\{ \mathbf{x} \, \Big| \sum\_{u \in G\_{\mathbf{r}}} so\_{uv} = \max \left( \sum\_{u \in G\_{i}} so\_{uv} \right) , \, 1 \le i \le R \right\}. \text{ A tie is broken as} \\ \text{choosing the error with the lowest number of machines.} \end{aligned}$$

by choosing the group with the lowest number of machines.

#### *3.4. GT Efficiency*

This work uses the GT efficiency measurement equation defined by Chandrasekharan and Rajagopalan [65], while letting *q* be equal to 0.5, as follows:

$$\eta = q\eta\_1 + (1 - q)\eta\_2 = q\frac{\varepsilon\_b}{\sum \quad Q\_i P\_i} + (1 - q)\left(1 - \frac{\varepsilon\_0}{\mu v - \sum \quad Q\_i P\_i}\right) \tag{9}$$

There are two concerns: the "1" in the group, and the "1" not in the group. The higher the GT efficiency, the more 1s are in the groups, and the fewer 1s are not in the groups. Thus, *eb* <sup>∑</sup> *NSiNOi* is the ratio of the number of 1s in the groups to the number of 1s and 0s in all groups, while *<sup>e</sup>*<sup>0</sup> *uv*−<sup>∑</sup> *NSiNOi* is the ratio of the number of 1s not in the groups to the number of 1s and 0s not in the groups.

#### **4. Numerical Examples**

*4.1. Example 1: Number of Operators at Each Workstation Decided*

On an assembly line, 500 operators must be assigned to workstations (Figure 6).

**Figure 6.** Assembly line with seven workstations.

Let

and

$$\begin{aligned} A &= \begin{bmatrix} a\_{11} & a\_{12} & a\_{13} \\ & a\_{22} & a\_{23} & a\_{24} & a\_{25} \end{bmatrix}, \\\\ \{t\_{ij}\} &= \begin{bmatrix} 10 & 30 & 24 & 0 & 0 \\ 0 & 5 & 40 & 20 & 15 \end{bmatrix} \end{aligned}$$

Let the work time per day be 8 h.

Thus,

$$
\begin{Bmatrix} \gamma\_{\bar{i}\bar{j}} \end{Bmatrix} = \begin{bmatrix} 48 & 16 & 20 & 0 & 0 \\ 0 & 96 & 12 & 24 & 32 \end{bmatrix}
$$

Thus, from Equation (6), *λline* = 5000/3, and from Equation (8),

$$
\begin{array}{cccc}
\{p\_{i\bar{j}}\} & ^{chop}\underline{\hspace{0.5cm}} & ^{c} & ^{104}\_{1} & \textbf{83} \\
\underline{\hspace{0.5cm}} & ^{c} & 17 & 138 & 69 & 52 \\
\end{array}
$$

The approximate operators needed at each workstation are as follows (Table 4).

**Table 4.** Approximate operators needed at each workstation.


The traditional step-by-step method is adjusted to obtain Table 5. Stations *a*<sup>11</sup> and *a*<sup>22</sup> are the bottleneck stations. Thus, if only one worker is assigned (the 498th one) to station *a*<sup>11</sup> or station *a*22, the assembly line cannot function. The production rate is increased if another worker is assigned to station *a*<sup>11</sup> and another is assigned to station *a*22. The final ALB solution is as follows (Table 5).

**Table 5.** Final ALB solution proposed by this study.


#### *4.2. Example 2: Auditing Production Rate*

The proposed method can be used to audit the production rate, which could be increased by rearranging the operators at each station. Consider this numerical example: An assembly line has seven workstations and 1100 workers (Figure 7).

**Figure 7.** Assembly line with seven workstations.

Let

$$A = \begin{bmatrix} a\_{11} & a\_{12} \\ a\_{21} & a\_{22} & a\_{23} \\ & & a\_{33} & a\_{34} \end{bmatrix}$$

and

$$\begin{aligned} \begin{Bmatrix} t\_{ij} \end{Bmatrix} &= \begin{bmatrix} 10 & 30 & \\ 20 & 40 & 15 \\ & 24 & 48 \end{bmatrix} \\\\ \begin{Bmatrix} p\_{ij} \end{Bmatrix} &= \begin{bmatrix} 70 & 180 & & \\ 120 & 200 & 80 & \\ & 150 & 200 \end{bmatrix} \end{aligned} $$

150 300

Let the work hours per day be 8 h. Thus, the following matrix is obtained:


From Equation (2), the following matrix is obtained:

 

 *λij* = ⎡ ⎣ 3360 2880 2880 2400 2560 3000 3000 ⎤ ⎦

where *λline* = 2400. If the workers are rearranged, the output rate is increased without any increase in cost. Equation (6) is used in the method proposed in this study: *λline* = 528,000/187. Thus, the following matrix is obtained:

ଵ ଽ ଼ ହ ସ ଷ ଶ ଵ ଵ ଵ ଵ ଶ ଵ ଷ ଵ ସ ଵ ହ ଵ ଵ ଵ ଼ ଵ ଽ ଶ ଶ ଵ ଶ ଶ ଶ ଷ ଶ ସ ଶ ହ ଶ ଶ ଶ ଼ ଶ ଽ ଷ ଷ ଵ ଷ ଶ ଷ ଷ ଷ ସ ଷ ହ


The traditional step-by-step method is adjusted to derive the final arrangement (Table 6):

**Table 6.** Final ALB solution proposed by this study.


#### *4.3. Example 3: Comparison with Boe and Cheng's Algorithm*

This example compares the proposed algorithm with Boe and Cheng's close neighbor algorithm [66]. The original part–machine incidence matrix is taken from Boe and Cheng as follows:

4.3.1. Machine Sorting

First, the relationship between the two machines must be obtained. For example, Equation (8) is used to obtain *b*37.

ଵ ଽ ଼ ହ ସ ଷ ଶ ଵ ଵ ଵ ଵ ଶ ଵ ଷ ଵ ସ ଵ ହ ଵ ଵ ଵ ଼ ଵ ଽ ଶ ଶ ଵ ଶ ଶ ଶ ଷ ଶ ସ ଶ ହ ଶ ଶ ଶ ଼ ଶ ଽ ଷ ଷ ଵ ଷ ଶ ଷ ଷ ଷ ସ ଷ ହ

݉ଷ ݉ ͳͳͳ ͳ ͳͳ ͳ ͳ ͳ ͳ ͳ ͳ ͳ ͳ ͳ ͳ ͳ ͳ ͳ ͳͳͳͳ <sup>൨</sup>

*b*<sup>37</sup> = 1 × 1 + 0 × 0 + 1 × 1 + ... + 0 × 0 = 6

Thus, the closeness matrix of skills is as follows.


In step 1 of machine sorting, the first machine, belonging to the first group, is *m*2. Thus, *G*<sup>1</sup> = {*m*2}. Ignoring *m*2, the next step is to form a new closeness matrix of skill *B u* = {*bkl*}, as follows:


The next machine is *m*14, which belongs to *G*<sup>1</sup> because *bm*14*m*<sup>2</sup> is the largest number in *B u*. Thus, *G*<sup>1</sup> = {*m*2, *m*14}. The new closeness matrix of skill *B u* = {*bkl*} is as follows:


The third machine is *m*1, which belongs to a new group, *G*2, because *bm*1*m*<sup>7</sup> is the largest number in *B u*, but machine *m*<sup>7</sup> does not belong to any group. Thus, *G*<sup>2</sup> = {*m*1}. Applying steps 2–3 in Section 3.3.1 gives four groups: *G*<sup>1</sup> = {*m*2, *m*14, *m*4, *m*18, *m*<sup>13</sup> }, *G*<sup>2</sup> = {*m*1, *m*7, *m*3, *m*8, *m*17, *m*5, *m*9},*G*<sup>3</sup> = {*m*15, *m*19, *m*11, *m*12, *m*<sup>16</sup> }, and *G*<sup>4</sup> = {*m*6, *m*10, *m*20}.

#### 4.3.2. Part Sorting

By using the operator sorting step in Section 3.3.2, the parts can be sorted. For example, *p*24, the ∑ *u*∈*G*<sup>1</sup> *soup*<sup>24</sup> = 5, and 2, 1, and 0 in *G*2, *G*3, and *G*4, respectively. Thus, *p*<sup>24</sup> belongs to *G*1.

Finally, the four groups are ordered (Table 7).

**Table 7.** The four groups in example 3.


The sorted incidence matrix is as follows:

ଵଵ ଵ ଵ ଶ ଶ ଷ ଵ ହ ଷ ଵ ଵ ଶ ଶ ଶ ଶ ଷ ଷ ଵ ଽ ସ ଶ ଶ ଷ ଷ ଷ ଵ ଼ ଵ ଵ ଶ ଶ

Boe and Cheng's solution of this example is as follows:

ଷ ଷ ଶ ଷ ସ ଷ ହ ଵ ଽ ସ ଵ ଶ ଵ ଶ ଼ ଷ ଷ ଵ ହ ଵ ହ ଵ ଶ ଶ ହ ଶ ଷ ଽ ଷ ଵ ଵ ଶ ଵ ଶ ଵ ଷ ଵ ଼ ଶ ସ ଶ ଵ ଼ ସ ଵ ଽ ଶ ଷ ଶ ଵ ଶ ଶ

The GT efficiency of the proposed algorithm is obtained using Equation (9), as follows:

$$\eta = 0.5 \times \frac{33 + 40 + 31 + 14}{5 \times 9 + 7 \times 11 + 5 \times 9 + 3 \times 6} + 0.5 \times \left( 1 - \frac{35}{20 \times 35 - \left( 5 \times 9 + 7 \times 11 + 5 \times 9 + 3 \times 6 \right)} \right) = 78.5\%$$

In the same way, GT efficiency is calculated using Boe and Cheng's algorithm, as follows:

$$\eta = 0.5 \times \frac{40 + 18 + 33 + 22}{7 \times 11 + 3 \times 8 + 5 \times 9 + 5 \times 7} + 0.5 \times \left(1 - \frac{40}{20 \times 35 - \left(7 \times 11 + 3 \times 8 + 5 \times 9 + 5 \times 7\right)}\right) = 77.4\%$$

#### *4.4. Example 4: Reducing Training Time*

This example demonstrates how to find an operator with the training and skills that matches a workstation's skill requirements so that the operator does not need to take any training courses. Consider the problem displayed in Figure 3 (Section 2). The proposed algorithm is used to obtain a sorted skill-operator matrix, as follows:

A manager could use the sorted skill–operator matrix as a database for choosing an operator with the right training and skills. Based on the database, the manager could quickly rotate a suitable worker to a workstation. Based on the sorted skill–operator matrix *SOsort*, the manager could choose operator *o*<sup>11</sup> for the workstation that requires skill *s*1; *o*2, *o*12, and *o*<sup>14</sup> because of *s*2; *o*4, *o*7, *o*13, and *o*3, because of *s*3; *o*<sup>1</sup> and *o*<sup>6</sup> because of *s*4; *o*<sup>15</sup> and *o*<sup>9</sup> because of *s*5; *o*<sup>5</sup> and *o*<sup>8</sup> because of *s*6; and *o*<sup>10</sup> because of *s*7, as follows:

Operators in the same skills group can support other operators in the group.

#### **5. Results and Discussion**

In example 1, stations *a*<sup>23</sup> and *a*<sup>24</sup> are the bottleneck stations. However, only one operator, namely the 500th one, can be assigned. Assigning the 500th operator to either station *a*<sup>23</sup> or *a*<sup>24</sup> cannot increase the output of the assembly line. This study suggests that the 500th worker should be laid off or shifted to another department. Alternatively, a 501st operator could be hired, which would allow another worker to be assigned to station *a*<sup>23</sup> and another to be assigned to station *a*24. This example demonstrates the proposed method in detail. The maximum output rate and the number of workers at each workstation can be decided. It also shows that multiple workstations can simultaneously become bottleneck stations.

In example 2, stations *a*<sup>12</sup> and *a*<sup>23</sup> are the bottleneck stations, and only one operator can be assigned. This study suggests that the 1100th worker be laid off or shifted to another department. Alternatively, a 1101st operator could be hired, which would allow another worker to be assigned to station *a*<sup>12</sup> and another to be assigned to station *a*23. With these solutions, a higher production rate than in the original arrangement (2816 > 2400) is achieved without any increase in cost. It shows that a manager can improve productivity by rearranging the number of workers at each workstation, while not increasing costs. It will also cut down on worker complaints because the workers' workloads are balanced.

In example 3, the proposed algorithm is used to categorize four groups, as shown in Section 4.3. Boe and Cheng's algorithm is compared with the GT efficiency measure proposed by Chandrasekharan and Rajagopalan. The GT efficiency of the proposed algorithm is 78.5%, higher than the 77.4% of Boe and Cheng's algorithm. Thus, the proposed GT algorithm works and can perform better than other algorithms in the literature. Workers are assigned to the workstation that best uses their skills, so they feel that their contribution is meaningful and more substantial. Workers at the same workstation also share more skills in common. This makes it easier for the manager to alternate workers.

In example 4, operators within the same skill group can support each other, giving managers more worker rotation solutions. With workers suggested by the sorted skill– operator matrix, managers can quickly match suitable operators with workstations so that productivity remains high. Using the matrix to match workers' training and skills to workstations' skill requirements also enables managers to reduce operator training time and cover workers' absences. For example, if *o*<sup>4</sup> is absent or busy, then *o*<sup>15</sup> can support or share the workload without needing too much training. The risk of an unexpected worker absence or resignation can be mitigated by alternating the workers according to the sorted skill–operator matrix.

#### **6. Conclusions**

Company managers aim to achieve sustainable operations, part of which is managing labor costs. To reduce labor costs, an assembly line manager increases the line's efficiency and reduces the training time needed by operators. Inadequate training or distractions,

such as an operator focused on a personal problem, can be dangerous, something a plant manager seeks to eliminate, while worker absences hurt efficiency and increase costs. This paper presents an algorithm that helps managers keep assembly lines balanced by arranging operators, so a line's workstations achieve the same output rates. The first example demonstrates how the proposed method can decide the maximum output rate and the number of workers at each workstation given a constant number of operators, especially when setting up a new plant. The second example demonstrates how the method can audit the existing number of operators at each workstation. Using the method to readjust the number of operators at a workstation can increase the output rate, without raising costs. As indicated in the numerical examples, another challenge is that multiple workstations can simultaneously become bottleneck stations. Knowing the number of operators at each workstation, the proposed cluster algorithm can group operators with similar skills and training, which reduces the training time. The third example demonstrates how the proposed clustering algorithm performs better than a previously described algorithm. Within the same group, operators from one workstation can also support operators at other workstations, as demonstrated in example 4. The clustering algorithm could be used by the human resources department to make a database that includes all of the company employees. When an employee is absent, a manger could quickly find a replacement, mitigating the risk of an unexpected absence. The manager could also use such a database to determine which employees should take a training course. By not having all of the employees take the training course, training costs could be reduced. Future applications of the proposed method could include the creation of a course-faculty database, where it could be used to find a replacement instructor for a class when the assigned instructor is absent.

Traditional GT considers only machines and parts, but this study uses the concept of GT to reduce the training time of workstation operators so that managers can achieve sustainable operations. The sorted skill–operator matrix in this study can, similar to a database, provide a manager with suggestions. However, most managers prefer a direct solution to a problem; they want the specific worker whose training and skills match the skill requirements of the problem workstation. Thus, future studies could develop an algorithm to assign specific workers to specific workstations by using a sorted skill– operator matrix, especially for large factories with multiple assembly lines.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Nomenclature**



#### **References**


### *Article* **A Novel Methodology for Simultaneous Minimization of Manufacturing Objectives in Tolerance Allocation of Complex Assembly**

**Lenin Nagarajan 1,\*, Siva Kumar Mahalingam 1, Sachin Salunkhe 1, Emad Abouel Nasr 2, Jõao Paulo Davim <sup>3</sup> and Hussein M. A. Hussein 4,5**


**Abstract:** Tolerance cost and machining time play crucial roles while performing tolerance allocation in complex assemblies. The aim of the proposed work is to minimize the above-said manufacturing objectives for allocating optimum tolerance to the components of complex assemblies, by considering the proper process and machine selections from the given alternatives. A novel methodology that provides a two-step solution is developed for this work. First, a heuristic approach is applied to determine the best machine for each process, and then a combined whale optimization algorithm with a univariate search method is used to allocate optimum tolerances with the best process selection for each sub-stage/operation. The efficiency of the proposed novel methodology is validated by solving two typical tolerance allocation problems of complex assemblies: a wheel mounting assembly and a knuckle joint assembly. Compared with previous approaches, the proposed methodology showed a considerable reduction in tolerance cost and machining time in relatively less computation time.

**Keywords:** tolerance allocation; machine and process selection; heuristic approach; univariate search method; whale optimization algorithm

#### **1. Introduction**

All aspects of manufacturing such as machine investment cost, manufacturing cost, the functionality of the product, quality of manufacturing, and the reliability of the product directly connect with tolerance allocation. Hence, most of the researchers are still focusing on this research topic to improve manufacturing efficiency. It involves the allocation of critical dimensions of an assembly, known from the product's functional requirements. As per the literature, more equations are available based on the combinations of the component's tolerance values; however, few equations provide better results. Tolerance allocation requires discovering the best possible combination of the component's tolerances by considering the manufacturing objective(s) and the associated constraints. The researchers proposed different tolerance allocation strategies in different periods. The summary of those strategies is explained here. Further, the detailed comparison of the strategies is given in the Supplementary File as Table S1.

**Citation:** Nagarajan, L.; Mahalingam, S.K.; Salunkhe, S.; Nasr, E.A.; Davim, J.P.; Hussein, H.M.A. A Novel Methodology for Simultaneous Minimization of Manufacturing Objectives in Tolerance Allocation of Complex Assembly. *Appl. Sci.* **2021**, *11*, 9164. https://doi.org/10.3390/ app11199164

Academic Editor: Vladimir Modrak

Received: 6 September 2021 Accepted: 27 September 2021 Published: 2 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Meta-Heuristic Method (MHM): In this method, the tolerance accumulation model, namely the root sum square (RSS) method [1–6], was used to calculate the assembly tolerance value by summing the randomly picked discrete tolerance values. Initially, the discrete tolerance values were identified by splitting the process tolerance limits. Most of the researchers had applied the MHM such as simulated annealing algorithm [7] and genetic algorithm [8–17] for this purpose.

Heuristic Method (HM): The application of this method is seldom found in literature, since the usage of thumb rules, previous experiences, and standards [18,19] in solving the tolerance allocation problems. Nevertheless, the alternative methods, namely the branch and bound algorithm [20] and design of experiments [21], were used to identify the effective tolerance allocation models. Further, a new method was developed by integrating HM with Tabu search for optimal tolerance allocation and subsequent manufacturing cost reduction [22]. Armillotta (2020) [23] minimized the manufacturing cost of the mechanical assemblies by properly allocating the tolerances and choosing the right dimensional properties. Korbi et al. (2021) [24] proposed a computer-aided design model for analyzing the tolerance in manufacturing the mechanical assemblies.

Discrete Cost Function (DCF) and Continuous Cost Function (CCF) Models: The researchers have used different cost function models [25–30] in the various periods to evaluate the manufacturing cost. These models were classified as DCF and CCF based on their nature. However, the researchers mainly used the CCF model, since it yields closed-form solutions to the tolerance allocation problems. On the other hand, the DCF models [4,26,31] were not preferred due to the model fitting errors during the manual formulation. Further, several studies have been carried out by considering the objective function as the sum of quality loss (as per Taguchi quality loss concept) and manufacturing cost [1,8,32–41].

Simple, Complex, and Non-Linear Assembly: The researchers used different tolerance allocation methods to obtain a solution concerning product type. For instance, the LMM was only used for simple products [1,8,42,43], which have only two mating components. Several works have been reported on finding the optimum allocated tolerance values for the components of complex assembly [44–51]. A limited number of authors have concentrated on non-linear assembly products that consist of more than two components [20,26,38].

Alternative Process Selection (APS): In practice, it is possible to produce components using more than one alternative process. It is necessary to select the proper process for the correct component to reduce the manufacturing cost. Some authors have considered alternative process selection (i.e., every combination of the process has a feasible tolerance range, and for a given process combination, the cost of machining is the function of the tolerance value) for optimum tolerance allocation of both simple and complex assemblies [4,10,13,15]. Kumar et al. (2009) [52] dealt with the distribution of tolerance on the component dimension of a complex assembly with alternative process selection. Authors have attempted to reduce the manufacturing cost of a product using the Lagrange multiplier method for complex assemblies with the bottom curve follower method.

Alternative Machine Selection (AMS): It is possible to reduce the manufacturing cost of a product by choosing the suitable machine for the correct process. Several researchers have considered machine selection as one of the criteria for minimizing manufacturing costs in recent years. However, very few authors [9,17,34] have discussed alternative machine selection for optimum tolerance allocation, and even fewer studies [34] on minimizing both cost and machining time with process and machine selection for optimum allocation of tolerance have been reported in the literature.

From this literature review, it appears that no significant effort has been made to consider machine time as an objective in optimum tolerance allocation, even though machine time is a crucial manufacturing parameter. Therefore, in the present study, both tolerance cost and machining time are optimized. Realizing the complexity of the problem, a combined heuristic and univariate search method is introduced to select the best machine

for each process and the optimum tolerance for each component using alternative process and machine selection, with available machine time as a constraint.

#### **2. Problem Definition**

The biggest challenge facing today's manufacturing companies is to reduce production costs while maintaining better quality and higher productivity. Tolerance allocation plays a vital role in achieving those goals. The production of a component involves selecting processes as well as selecting machines. These decisions directly influence the allocated tolerance values of components. The optimum allocated tolerance values govern the manufacturing costs and machining time of the product. Therefore, an operation may be possible with multiple alternatives, and as such, is treated as a non-polynomial hard problem.

#### **3. Mathematical Formulation**

The proposed work aims to simultaneously minimize tolerance cost and machining time, represented in Equation (1). Tolerance cost (*TCi*) and machining time (*MTik*) are calculated using Equations (2) and (3). Equations (4) and (5), respectively, determine the critical dimension of the sub-assembly (*Y*) and its tolerance (*tY*). Machine engagement time (met) is estimated using Equation (6). The constraints considered in this work are expressed by Equations (7)–(9). The allocated tolerance (*ti*) should be within the process limits, represented by Equation (7). The calculated sub-assembly tolerance within the given sub-assembly tolerance; this constraint is given by Equation (8). The total individual machine engagement time to manufacture the product should be less than the given available machine time (amt) represented in Equation (9). In addition, it is assumed that the following data are known well in advance before manufacturing the product:


$$Z = \min(TC, MT) \tag{1}$$

$$TC\_{ijk} = \eta\_{jk} \left( A\_j + \frac{B\_j}{t\_{ijk}} \right) \tag{2}$$

$$MT\_{ijk} = \eta\_{jk} \left( X1\_j + \frac{Y1\_j}{t\_{ijk}} \right) \tag{3}$$

$$Y = \sum\_{l=1}^{nc} \pm d\_l \tag{4}$$

$$t\_Y = \sum\_{i=1}^{no} t\_{ijk} \tag{5}$$

$$met\_k = \sum\_{i=1}^{no} MT\_{ijk} \tag{6}$$

$$t\_{\text{min}} \stackrel{\ast}{\leq} t\_i \stackrel{\ast}{\leq} t\_{\text{max}} \tag{7}$$

$$t\_a \underline{\chi} \ge t\_Y \tag{8}$$

$$amt\_k \ge met\_k\tag{9}$$


#### **4. Methodology**

The proposed method consists of two stages: (i) selection of the best machine for each process by applying a heuristic approach; (ii) selection of the best process and optimum allocated tolerance for each component using combined whale optimization algorithm and univariate search method. In the first stage, the process tolerance is divided into nd number of discrete values using Equation (10) and the allocated tolerance (*tejk*) is calculated using Equation (11). The tolerance cost (*TCejk*) and machining time (*MTejk*) for *tejk* are calculated using Equations (2) and (3), respectively. Nagarajan et al. (2018) explained that the distance method is used to combine the two different objective functions into a single one. For each discrete value, points are plotted on a graph where the *x*-axis and *y*-axis represent tolerance cost (*TCejk*) and machining time (*MTejk*), respectively. Assuming point (*x*1, *y*1) as the origin and point (*x*2, *y*2) as discrete tolerance cost and machining time, and substituting (*x*1, *y*1) as (0,0) and (*x*2, *y*2) as (*TCejk*, *MTejk*) in Equation (12), then the distance equation becomes Equation (13). The detailed steps of the heuristic and univariate search methods are shown in Figures 1 and 2. The pseudocode for the combined whale optimization algorithm and a univariate search method is presented in Section 5.

$$td\_{\dot{j}} = \frac{t \text{max}\_{\dot{j}} - t \text{min}\_{\dot{j}}}{nd} \tag{10}$$

$$t\_{ejk} = t \text{min}\_{\rangle} + (e - 1)td\_{\rangle} \tag{11}$$

where *e* is the index for discrete point of tolerance and takes from 1,2,3 . . . *nd*.

$$\text{dis} = \sqrt{\left(\mathbf{x2} - \mathbf{x1}\right)^2 + \left(y2 - y\mathbf{1}\right)^2} \tag{12}$$

$$\text{dis}\_{\varepsilon \bar{\text{jk}}} = \sqrt{\left(T\mathbb{C}\_{\varepsilon \bar{\text{jk}}} - 0\right)^2 + \left(MT\_{\varepsilon \bar{\text{jk}}} - 0\right)^2} = \sqrt{T\mathbb{C}\_{\varepsilon \bar{\text{jk}}}^2 + MT\_{\varepsilon \bar{\text{jk}}}^2} \tag{13}$$

**Figure 1.** Heuristic approach to determine the best machine for each process.

**Figure 2.** Flow chart of univariate search method.

#### **5. Numerical Illustration**

The proposed method is initially implemented in the existing problem (wheel mounting assembly) discussed by Geetha et al. (2013) [34] to show the method's effectiveness in case study 1. Later, it is implemented in knuckle joint assembly in case study 2. Case study 1: Wheel Mounting Assembly (WMA)

The components of the wheel mounting assembly are given in the Supplemental File as Figure S1. The operations required to manufacture the components of the assembly are illustrated in Figure S1, starting from *O*1 to *O*8. The feasibility of performing the operations through the processes from P1 to P5 is given in Table S2. The infeasibility of performing the operation using the specific process is marked as 0 . The same way, the feasibility of using machines for the specific processes is also given in Table S2. Further, the cost and time function constants are represented in Table S3. The first stage of the proposed work, the procedure to determine the best machine for each process, is discussed below.

Figure 3 represents the graphical representation of possible alternative processes and machines to carry out all the sub-stage operations to complete the manufacturing of the product. In the first stage, the best machine is selected for each process by implementing the heuristic approach.

**Figure 3.** Possible processes and machines for each operation/sub-stage to manufacture WMA.

Using Equations (2), (3), (10), (11) and (13), discrete tolerance (*tdj*), tolerance (*tejk*), tolerance cost (*TCejk*), machining time (*MTejk*), and distance values (*disejk*) are calculated for each machine and are shown in Table 1. Figure 4 is constructed for a discrete tolerance value of 0.01 mm (*tejk* = 0.01 mm), assuming *TCejk* on the *x*-axis and *MTejk* on the *y* axis. The distance between the origin and the discrete point of the machine represents the critical factor that decides on selecting the best machine. The value of less distance is more likely to be the best machine for the process.

**Table 1.** Best machine for process number 1.


**Figure 4.** Selection of the best machine based on discrete tolerance.

*TCe*11—tolerance cost of *e*th discrete tolerance using process number 1 on machine number 1; *MTe*13—machining time of *e*th discrete tolerance using process number 1 on machine number 3.

As shown in Figure 4, it is clear that *dis*-*M*1 is less than *dis*-*M*3; therefore, machine 1 is considered the best machine for process 1 for achieving a discrete tolerance value of 0.01 mm. The sum of the difference between minimum discrete distance (*min*(*disijk*)) and discrete distance (*disijk*) to that particular machine is calculated to select the best machine suitable for all discrete tolerances of the process tolerance. As shown in Table 1, it can be concluded that machine number 1 is the best machine for process number 1, shown graphically in Figure S2 (Supplementary File). Similarly, the best machine for the other processes is calculated and shown in Figure S2. It is clearly understood that machine numbers 2, 1, 3, and 4 are the best machines for process numbers 2, 3, 4, and 5, respectively. After implementing the heuristic approach, alternative machines are removed for each

process; instead, the best machine is selected, which is shown in Figure 5.

In the second stage, since the same allocated tolerances reported in Geetha et al. (2013) have been used for demonstration purposes, only univariate search method is implemented to get the best process, minimum tolerance cost, and minimum machining time for each substage/operation of WMA. The univariate search method proposed in this work provides the best results in relatively less computation time than the existing method. The search space of the existing method proposed in Geetha et al. (2013) contains 403,200 combinations, whereas, in the proposed method, there are only 11 possible combinations in the search space. Equations (14) and (15) compute the possible combinations in the existing and proposed methods, respectively. Table 2 represents the 11 possible combinations obtained using the univariate search approach, and the first combination is shown as yellow shaded text in Figure 5.

$$\begin{array}{c} ncs\_{\mathfrak{c}} = \prod\_{i=1}^{n0} \left( \sum\_{j=1}^{np} nm\_{ij} \right) \\ = (2+2+2)(2+3+2)(2+2)(3+3)(2+2)(2+3)(3+2)(2+2) \\ = 6 \times 7 \times 4 \times 6 \times 4 \times 5 \times 5 \times 4 = 403,200 \end{array} \tag{14}$$

$$nns\_{\mu} = 1 + \sum\_{i=1}^{n\nu} np\_i - n\nu = 1 + 3 + 3 + 2 + 2 + 2 + 2 + 2 + 2 - 8 = 11\tag{15}$$

**Figure 5.** Best machine for process of each operation/sub-stage after implementation of the first stage.

Using Equations (2) and (3), the tolerance cost and machining time are calculated for each operation with alternative processes and its corresponding best machine is shown in Figure 5. In all 11 combinations of alternative processes, the same tolerance mentioned in case 1 and case 2 by Geetha et al. (2013) is assumed for each operation. The alternative process and its corresponding machine-allocated tolerance of each sub-stage/operation and its tolerance cost and machining time are presented in Tables 2 and 3. The best result is shown as shaded text in both Tables 2 and 3.

Figures 6 and 7 show a comparison of the tolerance cost and machining time for the existing problem presented in Geetha et al. (2013) and the proposed method; it can be seen that the proposed method works well, considering tolerance cost, machining time, and both as objective functions.


**Table 2.** Result of univariate search method for the existing problem.

**Figure 6.** Comparison of tolerance cost of the proposed method with the existing problem presented in Geetha et al. (2013) MC: tolerance cost; EM10 & EM11: existing method, as per data presented in Tables 10 and 11 of Geetha et al.; PM10 & PM11: proposed method, as per data presented in Tables 10 and 11 of Geetha et al. (2013).

*Appl. Sci.* **2021**, *11*, 9164


**Figure 7.** Comparison of machining time of the proposed method with the existing problem presented in Geetha et al. (2013).

Table 4 represents the % of savings in both tolerance cost and machining time as compared to the existing method and the proposed method.

**Table 4.** % Savings in tolerance cost and machining time.


Case Study 2: Knuckle Joint Assembly

In practice, the availability of machine time will restrict the selection of the machine for performing a process, which will influence the tolerance cost and machining time. Therefore, in this work, the available machine time is considered as a constraint in selecting the machine. The methodology is demonstrated using a knuckle joint assembly (Figure 8), consisting of six components performed in ten sub-stages. Two critical dimensions are considered for the proper functioning of the product. It is assumed that nine processes are performed using ten machines. Equations (16) and (17) are used to determine the critical dimension, and the tolerances of the critical dimensions are estimated using Equations (18)–(21). Table 5 shows the dimensions, sub-stages, tolerance symbols, and tolerance stake-up of the knuckle joint assembly.

$$Y1 = d5 - d1 - d2 - d3 - d4 \tag{16}$$

$$
\Upsilon 2 = d 2 - d 6 \tag{17}
$$

$$t\_{Y1} = t\_{d5} + t\_{d1} + t\_{d2} + t\_{d3} + t\_{d4} \tag{18}$$

$$t\_{Y1} = t\_{O7} + t\_{O8} + t\_{O1} + t\_{O2} + t\_{O3} + t\_{O4} + t\_{O5} + t\_{O6} \tag{19}$$

$$t\_{Y2} = t\_{d2} + t\_{d6} \tag{20}$$

$$t\_{Y2} = t\_{O2} + t\_{O3} + t\_{O9} + t\_{O10} \tag{21}$$

$$\mathbf{C} = \sum\_{i=1}^{n0} T \mathbf{C}\_i \tag{22}$$

$$T = \sum\_{i=1}^{no} MT\_i \tag{23}$$

**Figure 8.** Knuckle joint assembly (KJA).


Table 6 shows the cost and time function constants of each process of the knuckle joint assembly. The possible operation–process and process–machine feasibility matrix and the available time of the individual machines are presented in Tables 7 and 8, respectively.

**Table 6.** Cost and time function constants of each process of knuckle joint assembly.



**Table 7.** Possible operation—process and process—machine feasibility matrix of knuckle joint assembly.

**Table 8.** Available machine time of each machine in the knuckle joint assembly.


As per stage 1, the best machine for each process is determined using the heuristic method and is presented in Table 9 and Figure S3 (in Supplementary File).

**Table 9.** Best machine for each process of the knuckle joint assembly.


In the second stage, the whale optimization algorithm is implemented with a univariate search method. The optimum allocated tolerance for each sub-stage/operation, tolerance cost, machining time, and best process is obtained by assuming 100 whales, 100 iterations as stopping criteria with 20 runs. The parameters involved in the algorithm are given in Table 10. Further, the pseudocode of this algorithm is given in the Supplementary File.

**Table 10.** Terms involved in whale optimization algorithm.


The convergence plot for tolerance cost and machining time for stopping criteria considered as 100 iterations are shown in Figure 9. The Pareto optimal solution for a sample run is shown in Figure 10. The tolerance cost and machining time for 31 runs are presented in Table 11. Out of these 31 runs, the best value is calculated using EDAS and CODAS multi-criteria decision-making techniques implemented by Adali & Tu¸s (2019) [53] The appraisal score (APS) and the assessment score (AS) obtained by implementing EDAS and CODAS method are presented in Table 11. Run number (*R*.*No*.) 11 has the highest APS value in EDAS, and run number 7 has the highest AS value in CODAS (Highlighted in Table 11). The optimum allocated tolerances of each sub-stage are presented for both methods and are listed in Table 12, along with the tolerance cost and machining time for the required sub-assembly tolerances of 0.55 mm and 0.22 mm. The optimum total tolerance cost and total machining time are presented in Table 13 for the given sub-assembly tolerance values with the constraint of available machine time.

**Figure 9.** Convergence plot for tolerance cost and machining time.

**Figure 10.** Optimum Pareto front solutions.

Table 14 shows the machine engagement times of the individual machines to manufacture the knuckle joint assembly within the available machine time.

For supporting the proposed method, the statistical analysis for EDAS and CODAS methods are executed through Minitab software. The statistical analysis results and probability plots for both the methods are presented in Figure 11. The probability values are 0.329 and 0.231 for the EDAS and CODAS methods, respectively. Since the value of probability is greater than 0.005 in both cases, it is clearly understood that the results obtained in 31 runs are from normally distributed data.


**Table 11.** Appraisal (APS) and assessment (AS) scores of EDAS and CODAS.

**Table 12.** The optimum allocated tolerance of knuckle joint assembly.


*O*.*No*.—operation/sub-stage number; *P*.*No*.—process number; *M*.*No*.—machine number, *Tol*.—allocated tolerance; *TC*—tolerance cost; *MT*—machining time.


**Table 13.** Sub-assembly tolerance, total cost, and machining time.

*ta*1—Calculated sub-assembly 1 tolerance. *ta*2—Calculated sub-assembly 2 tolerance.

**Table 14.** Machine engagement time (*met*) of individual machine.


**Figure 11.** Statistical analysis and probability plot for EDAS and CODAS methods.

#### **6. Conclusions**

Most previous studies on tolerance allocation problems concentrated on minimizing manufacturing costs, quality loss, or combining the two. Machining time, a vital manufacturing objective, has barely been contemplated. In this paper, the machining time was considered along with manufacturing cost in optimum tolerance allocation of complex assemblies, representing a more realistic product development scenario. Alternative machine and process selections with available machine time make this problem cumbersome

and complicated. Therefore, a new methodology was developed that applies a heuristic approach and combines whale optimization algorithm with a univariate search method. The total manufacturing cost and machining time of 36.3 USD and 57.74 min reported in this paper for wheel mounting assembly is 2.73% and 15.35% less than the problem dealt with in case 1 by Geetha et al. (2013). Similarly in case 2, there was 6.22% and 20.32% of savings in the tolerance cost and machining time reported by implementing the proposed method. The results presented in this paper demonstrate that the proposed method can reduce tolerance cost and machining time simultaneously with less computation time. The proposed method is also suitable for solving two- and three-dimensional problems. As a further extension of this work, the operation sequence, machine sequence, or both may be considered with additional objectives such as total investment cost of machines, idle time of machines, idle cost of machines, and the number of machines required to manufacture the product.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/app11199164/s1, Figure S1: Wheel Mounting Assembly (WMA) [Geetha et al. (2013)], Figure S2: Selection of best machine for process number 1 to 5 to manufacture WMA, Figure S3: Selection of best machine for process number 1 to 9 to manufacture KJA, Table S1: Summary of literature survey, Table S2: Feasibility Matrix for WMA [Geetha et al. (2013)], Table S3: Cost and Time Function Constants for WMA [(Geetha et al. 2013)] [54–69].

**Author Contributions:** Conceptualization and methodology by L.N. and S.K.M.; writing—original draft preparation by L.N. and S.K.M.; writing—review and editing by S.S., E.A.N., J.P.D. and H.M.A.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors extend their appreciation to King Saud University for funding this work through Researchers Supporting Project number (RSP-2021/164), King Saud University, Riyadh, Saudi Arabia.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare that they have no conflict of interest.

#### **References**


### *Article* **Harmony Search Algorithm for Minimizing Assembly Variation in Non-linear Assembly**

**Siva Kumar Mahalingam 1, Lenin Nagarajan 1,\*, Sachin Salunkhe <sup>1</sup> Emad Abouel Nasr 2, Jõao Paulo Davim <sup>3</sup> and Hussein M. A. Hussein 4,5**


**Abstract:** The proposed work aims to acquire the maximum number of non-linear assemblies with closer assembly tolerance specifications by mating the different bins' components. Before that, the components are classified based on the range of tolerance values and grouped into different bins. Further, the manufacturing process of the components is selected from the given and known alternative processes. It is incredibly tedious to obtain the best combinations of bins and the best process together. Hence, a novel approach using the combination of the univariate search method and the harmony search algorithm is proposed in this work. Overrunning clutch assembly is taken as an example. The components of overrunning clutch assembly are manufactured with a wide tolerance value using the best process selected from the given alternatives by the univariate search method. Further, the manufactured components are grouped into three to nine bins. A combination of the best bins is obtained for the various assembly specifications by implementing the harmony search algorithm. The efficacy of the proposed method is demonstrated by showing 24.9% of cost-savings while making overrunning clutch assembly compared with the existing method. The efficacy of the proposed method is demonstrated by showing 24.9% of cost-savings while making overrunning clutch assembly compared with the existing method. The results show that the contribution of the proposed novel methodology is legitimate in solving selective assembly problems.

**Keywords:** selective assembly; overrunning clutch assembly; univariate search method; harmony search algorithm

### **1. Introduction**

Product quality is the focus of any manufacturing process. In general, two or more components are assembled to create an assembly. The quality of the assembly depends on the quality of the components, which affects the product's functionality. Tolerance plays an essential role in the component's quality, deciding the fit between the mating parts. The components manufactured with closer tolerance make the precise assembly more suitable for functional requirements. Selective assembly is one of the feasible methods for making precise assemblies with lower manufacturing costs. The complete elimination or reduction of secondary operations by forming wide-tolerance components is the reason for lower manufacturing costs. In selective assembly, the components are grouped as bins based on the uniform grouping method, the equal probability method, or the uniform tolerance method, for example. According to the best bin combinations, the precise assemblies

**Citation:** Mahalingam, S.K.; Nagarajan, L.; Salunkhe, S.; Nasr, E.A.; Davim, J.P.; Hussein, H.M.A. Harmony Search Algorithm for Minimizing Assembly Variation in Non-linear Assembly. *Appl. Sci.* **2021**, *11*, 9213. https://doi.org/10.3390/ app11199213

Academic Editor: Vladimir Modrak

Received: 6 September 2021 Accepted: 28 September 2021 Published: 3 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

are made by mating the components randomly selected from the corresponding bins. Most of the research on linear/radial assembly have mainly focused on minimizing the objectives, namely, clearance variation or surplus parts. Research works focusing on non-linear assembly are seldom found. Moreover, in the existing literature, the tolerance of the components is usually considered for obtaining the best bin combinations rather than the dimension of components. Since different combinations of bins are possible, this problem environment can be treated as an NP-hard problem. Selection of process for making components from the given alternative processes also plays a vital role in further reducing the product's manufacturing cost.

#### **2. Related Research**

The literature relevant to the proposed work has been categorized as selective assembly and the harmony search algorithm. The literature related to both topics is detailed below.

#### *2.1. Selective Assembly*

Kern (2003) [1] proposed an approach for selecting assembly by considering variations in the dimensional distributions. Further, the author developed closed-form equations for the different techniques of selective assembly. Mease et al. (2004) [2] introduced a method to classify the components into various classes based on the dimensional variations. The assembly was made by pairing the selected components from the different classes. Further, the absolute and squared error loss functions were studied to select the optimum method. Matsuura et al. (2008) [3] explored a method to minimize clearance variation when pairing the components with less variance for making an assembly. Kwon et al. (2009) [4] studied the effectiveness of optimal binning strategies based on the squared error loss function by considering similar normal distributions of the dimensions of two different components. Matsuura et al. (2011) [5] developed the optimal mean shift to reduce surplus components while making selective assembly. The equal width concept was used for grouping the components into different bins. Three shifted means method was used for the fabrication of components.

Yue et al. (2014) [6] proposed a genetic algorithm to minimize the variation in clearance while mating the components available in the selective groups to manufacture hole and shaft assemblies. Babu and Asha (2014) [7] introduced a customized artificial immune system algorithm for mating the components available in selective groups to minimize assembly clearance variation. Further, they applied Taguchi's loss function to identify the deviation from the mean. The count of selective groups was also one of the important concerns in the proposed work. Ju and Li (2014) [8] solved selective assembly problems using a two-stage decomposition procedure. A two-component assembly system with unreliable Bernoulli machines and limited inventories was considered. The effectiveness of the proposed method was analyzed analytically.

Xu et al. (2014) [9] applied a new method with the combination of discarding and binning theorems to minimize the components' variations while making the selective assembly of a hard disk drive. The utilization rate of the components was high with regard to the number of assemblies made using the newly proposed method. Lu and Fei (2015) [10] proposed a grouping method in selective assembly to increase the success rate of manufacturing the assemblies and reducing the surplus parts. A genetic algorithm with a 2D chromosome structure was used for this purpose. The effectiveness of the proposed method was studied based on solving three different cases. Manickam and De (2015) [11] developed a genetic algorithm to identify the better combination of groups containing the components to create the maximum number of assemblies with less total cost. The small number of components with wider tolerance was made available in the individual groups. This was the uniqueness of the presented work. The model was developed using MATLAB software.

Babu and Asha (2015) [12] evaluated the losses in making assemblies by applying the symmetrical interval-based Taguchi loss function. Further, the dual objectives, namely, minimum clearance variation and minimum losses in making assemblies, were considered

to identify the better combination of the parts' groups using the sheep flock heredity algorithm. Ju et al. (2016) [13] studied the performance of the selective assembly system using Bernoulli machine reliability models. For this study, the assemblies made through two components with different qualities were considered. The two-stage decomposition procedure was applied to evaluate the effectiveness of the proposed models analytically. Malaichamy et al. (2016A) [14] developed software capable of applying all kinds of selective assembly methods to identify suitable solutions with the least possible time. For instance, equal area and width methods were used to minimize variation and scrap in selective assembly. The random assembly method was inbuilt in the software to select the best method for minimizing the manufacturing cost of the assembly. Since the software was developed using an optional programming language, changes to the input data were extremely easy, and, in no time, the result could be viewed graphically.

Malaichamy et al. (2016B) [15] introduced a software approach to help engineers to visualize the tolerance cost curve for a given process or a good combination of processes in no time. The tolerance cost was calculated using the reciprocal power cost model. The best process could be selected in this work from the given alternative processes for the known tolerance value. A simulated annealing algorithm was implemented to compute the best bin number and its combinations with almost nil surplus parts for the given dimensional distribution of components of a shaft housing assembly. Compared with random assembly, a good number of surplus parts were reduced by implementing an equal width selective assembly method. Liu and Liu (2017) [16] discussed the remanufacturing of the engine through the selective assembly concept. Remanufacturing is a sustainable strategy in which the group numbers and range of components in each group are not kept constant. The assembly accuracy was greatly increased using the proposed strategy. Chu et al. (2018) [17] developed a method for manufacturing gear reducers using a novel strategy called GAbased selective assembly. The backlash of the gear reducer was the major concern in this work to verify the meeting of assembly requirements.

#### *2.2. Harmony Search Algorithm*

Geem et al. (2001) [18] described the harmony search optimization algorithm's features, robustness, simplicity, and search efficiency in solving engineering problems. Geem et al. (2002) [19] proposed a harmony search optimization algorithm to solve the problems associated with pipe network design. Lee and Geem (2004) [20] described a novel structural optimization algorithm using the harmony search method. The effectiveness of the proposed algorithm was verified by identifying the solutions to different truss problems. Lee and Geem (2005) [21] solved different engineering optimization problems that included minimizing mathematical functions and optimizing structural problems using a harmonic-search-based method. Flow demand and low head were considered objectives in this work. Wang et al. (2009) [22] analyzed the harmony search algorithm's effectiveness with the colonal selection algorithm. The colonal search method was used to increase the harmony memory members in the harmony search method.

From the literature survey, it is clear that research works have been very limited in the manufacturing of non-linear assemblies with a closer assembly tolerance specification. Further, the usage of the harmony search algorithm for solving selective assembly problems is seldom found. The problem environment is discussed in the next section.

#### **3. Problem Environment**

In selective assembly, parts are manufactured with wider tolerance, measured, and partitioned into groups, and assemblies are made by assembling the components within the random combination of groups. Manufacturing tolerance, assembly specification, and the number of groups are the three main factors that play important roles in controlling surplus parts that affect the product's manufacturing cost. After parts are manufactured, it is tedious to obtain the best combination of groups for different assembly specifications. It is laborious work to compute the numbers of closer assembly for each bin number from the manufactured components. Alternative process selection for making components also makes the problem highly complicated. The problems described above are challenging to process/design engineers and reduce implementation in real situations. Moreover, the components involved in making a non-linear assembly require the individual mating of dimensions for each assembly. Compared to making linear assembly, this required additional computation effort is tedious.

#### **4. Solution Methodology**

The proposed problem environment can be solved in three stages. In the first stage, the best process for manufacturing each component of an assembly is obtained from the alternative processes using the univariate search method. In the second stage, 1000 simulated dimensions of each component are generated using MATLAB for the different combinations of alternative processes obtained from the first stage. Further, the components are partitioned into different bin numbers according to the equal area method, which is one of the techniques used in the selective assembly method. In the last stage, the harmony search algorithm is implemented to obtain the best bin combinations. Then, the non-linear assemblies are made by mating the components according to the best bin combinations with almost nil surplus parts. To show the effectiveness of the proposed method, the manufacturing cost of assemblies made both from the random assembly [23] and selective assembly are compared. This three-stage approach is explained in a detailed way in the numerical illustration section.

#### **5. Numerical Illustration**

Overrunning clutch assembly (OCA), dealt with by Ganesan et al. (2001) [23], shown in Figure 1, is considered an example product to show the efficiency of the proposed method. It consists of a hub, four rollers, and a cage. The nominal dimensions and their allocated tolerance; the minimum, maximum process tolerances; and the tolerance cost function constants of alternative processes of each component are listed in Table 1. The critical dimension Y is determined using Equation (1), and the accepted value of Y is assumed to be 7.0124 ± 2◦. Equation (2) represents the reciprocal tolerance cost function to calculate the component's tolerance cost. The cost function constants are taken from Ganesan et al. (2001) [23]. The cost of assembly based on allocated and maximum process tolerances is computed using Equations (3) and (4), respectively, and listed in Table 2.

**Figure 1.** Overrunning clutch assembly.

$$Y = \cos^{-1}\left(\frac{X1 + \left(\frac{X2 + X3}{2}\right)}{X4 - \left(\frac{X2 + X3}{2}\right)}\right) \tag{1}$$

$$C\_i = A\_i + \frac{B\_i}{t\_i} \tag{2}$$

$$\mathbf{C}\_{d} = \sum\_{i=1}^{nc} A\_{i} + \frac{B\_{i}}{t\_{ai}} \tag{3}$$

$$C\_m = \sum\_{i=1}^{nc} A\_i + \frac{B\_i}{t\_{mi}} \tag{4}$$


**Table 1.** Manufacturing details of overrunning clutch assembly.


*C.No.*—component name; *P.No.*—process number.

**Table 2.** Manufacturing cost of overrunning clutch assembly for *tai*, *tmin*, and *tmax*.


#### *5.1. Stage I*

The selection of the best process for each component from the given alternative processes using the univariate search method is illustrated in Figure 2. It is understood from Figure 2 that the minimum manufacturing cost of \$18.37 can be achieved through making the components *X1*, *X2*, *X3*, and *X4* using process combinations (3213) *P3, P2, P1*, and *P3*, respectively, which is nearly 24.98% less compared with the \$24.49 reported in the existing method dealt by Ganesan et al. (2001). It is also understood that the process combinations of 2213 and 1213 for components *X1*, *X2*, *X3*, and *X4*, respectively, can yield 17.52% and 9.92% savings in manufacturing cost. However, in a real situation, the savings may vary slightly because of surplus parts present in the selective assembly method.

**Figure 2.** Univariate search method to select the best process for each component. *X1-P1* indicates that *X1* component is produced by the *P1* process; TC1113 = \$26.56 indicates that components *X1*, *X2*, *X3*, and *X4* are manufactured using processes *P1, P1, P1*, and *P3*, respectively, and the total cost to manufacture the same will be \$26.56.

#### *5.2. Stage II*

As discussed in Section 4, 1000 random values have been generated for each component according to the mean (*μi*) and standard deviation (*σi*) presented in Table 3. The dimensional distribution of 1000 components of *X1*, *X2*, *X3*, and *X4* was generated using the normrnd (*C.No.*) function in MATLAB.


**Table 3.** Mean and standard deviation of components for different process combinations.

*P.No.*—process number; *PC1*—1st process combinations for components *X1*, *X2*, *X3*, and *X4*; *μi*—mean dimension of the components; *tmi*—tolerance of the *i*th component; *σi*—standard deviation of the *i*th component.

#### *5.3. Stage III—Implementation of HSA*

The harmony search algorithm (HSA) is a meta-heuristic algorithm, and it works based on the identification of good harmony by musicians through a continuous improvisation process. The HSA has the following advantages: (i) quick convergence, (ii) easy to adapt, and (iii) the least computational time. Further, from the literature survey, it is observed that the HSA has outperformed in solving complex optimization problems. Hence, the HSA has been used in this work. Table 4 presents the different terms used in the HSA, its equivalent term in both optimization problems and the present work formulation, and its range of values and examples. The schematic diagram shown in Figure 3 illustrates the implementation of the HSA to obtain the best bin combinations. The technical terms and their meanings used in Figure 3 are presented in Table 5. For demonstration purposes, the components are partitioned into five bins. The step-by-step procedure is given below.

**Table 4.** Representation of variables in HSA.


**Figure 3.** Implementation of HSA.


**Table 5.** Technical terms and their meanings.

Step 1: A random combination of bins for each component for the size of harmony memory 10 is generated (listed in Table 6).


**Table 6.** Initial harmony memory.


*H.No.*—harmony numbers; *cbX1*, *cbX2*, *cbX3*, and *cbX4* are the combinations of bins for *X1*, *X2*, *X3*, and *X4* components.




**Table 8.** Best combination of bins for the first iteration.

Step 4: A random number less than 1 is generated for each component of each harmony number (presented in Table 9 as *RHMCR*).

Step 5: If this number is less than or equal to *HMCR*, then a random number between 1 and *HMS* (*rcb*) is generated, and others are assumed to be zero. This is presented in Table 9.

**Table 9.** *RHMCR* and *rcb* values.


*RHMCR*—a random number decides either to accept or not accept the selection of the existing tune/bin combination of components for improvisation; *rcb*—a random number generated between 1 and hms for each variable.

Step 6: The *cbX1, cbX2, cbX3,* and *cbX4* values corresponding to *rcbX1, rcbX2, rcbX,* and *rcbX4* are taken from Table 6 and listed in Table 10. If the value of *rcb* is zero, then the corresponding harmony's *cbx* value is considered.

**Table 10.** Harmony after HMCR.


Step 7: A random number (*RPAR*) less than 1 is generated for each value of *RHMCR*, which is not equal to zero and is less than the *HMCR* value (listed in Table 11).

Step 8: Two random numbers, r1 and r2, within bin numbers, are generated for each value of RPAR, which is not equal to zero (presented in Table 11).

Step 9: New harmony, i.e., *cbX1*, *cbX2*, *cbX3*, and *cbX4*, is obtained wherever the bin is located within their bin combinations, according to r1 and r2 (presented in Table 12).

Step 10: Then, NoAs are obtained by mating the components randomly corresponding to the bin combinations given in *cbX1*, *cbX2*, *cbX3*, and *cbX4* (listed in Table 12).

Step 11:The selection and selected harmonies for the next iteration are presented in Tables 13 and 14.

Step 12: The steps from 4 to 11 are repeated to the specified number of iterations.


**Table 11.** *RPAR, r1*, and *r2* values.


*RPAR*—a random number decides the pitch adjustment (interchanging bin numbers within the component); *r1* and *r2*—two random numbers within the bin number (*nb*) to interchange the bin number to generate a new bin combination for a component.



**Table 13.** Selection of harmonies for next iteration.


*NS*—not selected; *SHMNo.*—selected harmonies for next iteration.


**Table 14.** Selected harmonies for next iteration.

**E** *NoA* for ± 1°

**Figure 4.** *Cont*.

(**d**). *NoA* for ±2°

**Figure 4.** Iteration number vs. *NoA* for various bin numbers—the 3213 process combination (**a**) *NoA* for ±0.5◦, (**b**) *NoA* for ±1◦, (**c**). *NoA* for ±1.5◦ and (**d**) *NoA* for ±2◦.

**Figure 5.** *Cont*.

**Figure 5.** *NoA* for various bin numbers and assembly specifications for the 3213 process combination (**a**) *NoA* for ±0.5◦, (**b**) *NoA* for ±1◦, (**c**). *NoA* for ±1.5◦ and (**d**) *NoA* for ±2◦.

#### **6. Results and Discussion**

An attempt has been made to make assemblies by considering the number of bins/ partitions, up to 9. Figure 6 reveals that in the equal area method for making non-linear assemblies, it is possible to make 996 assemblies out of 1000 components by partitioning them into 3 bins for the assembly specification of ±2◦. It is also understood that while increasing the partition number, there may be a 5.8% (938 assemblies) drop in producing the number of assemblies for the same assembly specification. In the meantime, the number of assemblies is reduced for the same 1000 components of *X1*, *X2*, *X3*, and *X4* while reducing the assembly specification for the same partition number (equal to 3). The assemblies are reduced from 996 to 392 for the assembly specification of ±2◦ to ±0.5◦.

**Figure 6.** *NoA* vs. bin number for various assembly specifications for the 3213 process combination.

Figures 7–9 represent the number of assemblies produced for various assembly specifications while changing the partition number for the same set of 1000 components produced based on process combination 3213. A similarity can be observed from the above figures. Except in three bin partitions, all other partitions of components and matching components according to the best bin combinations obtained through the HSA produced almost a very close number of assemblies for various assembly specifications. Figure 10 indicates that maximum assemblies are produced for various assembly specifications and process combinations while partitioning the components into three bins. The assembly specification value after ±1◦ in all the process combinations could produce an almost equal number of assemblies for bin number three. Table 15 represents the best bin combinations and their maximum number of assemblies for various process combinations.

**Figure 7.** Assembly specification vs. *NoA* for various bin numbers—the 3213 process combination.

**Figure 8.** Assembly specification vs. *NoA* for various bin numbers—the 1213 process combination.

**Figure 9.** Assembly specification vs. *NoA* for various bin numbers—the 2213 process combination.

**Figure 10.** Maximum *NoA* produced for different assembly specifications and various process combinations.


**Table 15.** Best bin combinations for various assembly specifications and process combinations.

*AS*—assembly specification.

#### **7. Conclusions**

This paper addresses a novel methodology by combining the univariate search method and the harmony search algorithm in selective assembly for making non-linear assemblies for various assembly specifications. The best processes for the different components of the assembly are selected from the known alternative processes using the univariate search method, and these components are grouped into 3 to 9 bins. Further, the best bin combinations for making assemblies to reduce manufacturing cost are obtained through the harmony search algorithm. In this work, the component's dimensions are directly considered for making assemblies from the best bin combinations rather than considering tolerances, as in the existing method. The proposed method is demonstrated on a non-linear overrunning clutch assembly and has proved its efficiency by saving 24.9% of manufacturing cost compared with the existing method for the best process combination of 3213.

**Author Contributions:** Conceptualization and Methodology by L.N. and S.K.M.; Writing—Original Draft Preparation by L.N. and S.K.M.; Writing—Review and Editing by S.S., E.A.N., J.P.D. and H.M.A.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors extend their appreciation to King Saud University for funding this work through Researchers Supporting Project number (RSP-2021/164), King Saud University, Riyadh, Saudi Arabia.

**Institutional Review Board Statement:** Not Applicable.

**Informed Consent Statement:** Not Applicable.

**Data Availability Statement:** Not Applicable.

**Conflicts of Interest:** The authors declare that they have no conflict of interest.

**Availability of Data and Material:** Not applicable.

#### **References**


### *Article* **A Hybrid Bat Algorithm for Solving the Three-Stage Distributed Assembly Permutation Flowshop Scheduling Problem**

**Jianguo Zheng and Yilin Wang \***

Glorious Sun School of Business and Management, Donghua University, Shanghai 200051, China; zjg@dhu.edu.cn **\*** Correspondence: yilinwang0319@163.com

**Abstract:** In this paper, a hybrid bat optimization algorithm based on variable neighbourhood structure and two learning strategies is proposed to solve a three-stage distributed assembly permutation flowshop scheduling problem to minimize the makespan. The algorithm is firstly designed to increase the population diversity by classifying the populations, which solves the difficult trade-off between convergence and diversity of the bat algorithm. Secondly, a selection mechanism is used to update the bat's velocity and location, solving the difficulty of the algorithm to trade-off exploration and mining capacity. Finally, the Gaussian learning strategy and elite learning strategy assist the whole population to jump out of the local optimal frontier. The simulation results demonstrate that the algorithm proposed in this paper can well solve the DAPFSP. In addition, compared with other metaheuristic algorithms, IHBA has better performance and gives full play to its advantage of finding optimal solutions.

**Keywords:** hybrid bat algorithm; optimization problem; the distributed assembly permutation flowshop scheduling problem; variable neighborhood descent

#### **1. Introduction**

The distributed assembly scheduling problem is a derivative and extension of the assembly scheduling problems. It is mainly used to produce multiple products assembled from different parts. In a decentralized and globalized economy, the current economic trend of customization and intelligent manufacturing has led to the rapid development of assembly production. Internationally, international companies with multiple production centers or plants are even more common. This shows that distributed and intelligent production plays a pivotal and irreplaceable role in the modern manufacturing industry. Currently, this type of scheduling is widely used in the supply chain and manufacturing industry. In order to maintain a competitive position in a rapidly changing market, managers must quickly make the right decisions about how to allocate work to plants and how to schedule each plant efficiently. Therefore, the distributed scheduling problem has become a hot research topic among scholars and researchers.

The distributed assembly permutation flowshop scheduling problem (DAPFSP) is of great importance to both the industrial industry and the research community. As we all know, DAPFSP consists of two phases, including production and assembly. In detail, the production phase corresponds to a distributed displacement flowshop scheduling problem, where *n* jobs are assigned to be processed in *f* plants with *m* identical machines in a flowshop line. Each job belongs to a certain product. The processes contained in each product cannot be interrupted or inserted by processes of other products during the production process. The assembly phase corresponds to the assembly flow shop scheduling problem, where *s* products are made in an assembly plant with a single assembly machine, and each product must be assembled after all processes produced in the production phase are completed.

Total completion time has been identified as a more relevant and meaningful goal in today's dynamic manufacturing environment when a batch of work needs to be completed

**Citation:** Zheng, J.; Wang, Y. A Hybrid Bat Algorithm for Solving the Three-Stage Distributed Assembly Permutation Flowshop Scheduling Problem. *Appl. Sci.* **2021**, *11*, 10102. https://doi.org/10.3390/ app112110102

Academic Editors: Vladimir Modrak and Zuzana Soltysova

Received: 17 September 2021 Accepted: 22 October 2021 Published: 28 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

as quickly as possible. This allows the objective function to minimize the flow of work and efficiently improve resource utilization. DPFSP is an NP-hard with *m* greater than or equal to 2, and is a derivation and extension of the traditional permutation flow shop scheduling problem with this criterion. Similarly, DAPFSP can be seen as an extension of DPFSP, which is more complex than DPFSP. Therefore, the DAPFSP with completion time as the objective function is also a strong NP-hard problem. In recent years, as a simple yet efficient method, the bat algorithm has solved many scheduling problems in academia, including the job shop scheduling problem [1], the open shop scheduling problem [2], the task scheduling problem [3], and the production scheduling problem [4], by continuously performing local search in the neighborhood of the current solution to obtain the optimal solution.

In this paper, an improved hybrid bat algorithm, IHBA, is proposed for solving the DAPFSP with maximum makespan based on the original bat algorithm. The main contributions of this paper are:


This paper is organized as described below. Section 2 reviews the closely related literature. Section 3 presents the DAPFSP problem and its mathematical model. In Section 4, the hybrid bat algorithm proposed in this paper is described in detail. In addition, the parameter calibration and simulation calculation results are analyzed in Section 5. In the last section, the paper is summarized, and future work is provided.

#### **2. Literature Review**

It is well known that the distributed permutation flowshop scheduling problem (PFSP) is one of the most well-known production scheduling problems. It has a strong engineering background by having the same work order on all machines. The problem has been shown to be a strong NP-hard problem when more than two machines are involved.The problem of distributed assembly permutation flowshop scheduling problem (DAPFSP) has been gaining attention and popularity among scholars and researchers. The most leading exponent is Hatami et al. [5], whose work is very enlightening. He was the first to optimize and improve the DAPFSP for modeling and studying complex supply chains in 2013. They also introduced the mixed integer linear programming model, proposed three construction algorithms, tested the variable neighborhood descent (VND) algorithm, and demonstrated the good performance of the VND algorithm in solving this scheduling problem. Similarly, he went on to expand on his previous paper by adding time for sequence-related setups in both the production and assembly stages [6]. In 2015, Hatami et al. [7] considered a distributed assembly scheduling problem with sequence-dependent setup times and the goal of minimizing production time. The setup times of both phases are sequence-dependent, which also lays the theoretical foundation for the sequences in this paper. Heuristics and meta-heuristics are also proposed to solve it. Based on his literature and views, many of these ideas and opinions have been studied and explored in greater depth by many scholars. Based on these rather cutting-edge ideas, many scholars have built on and deepened their research, and Ying et al. [8] extended the DAPFSP for flexible assembly and sequence-independent setup times in supply chains. Using completion time as an optimization criterion, construction heuristic and custom meta-heuristic algorithms are proposed to solve this emerging scheduling extension. Gonzalez-Neira et al. [9] studied a stochastic version of the DAPFSP and proposed a hybrid algorithm to solve the stochastic problem, integrating biased randomization and simulation techniques. Wang [10] introduced a just-in-time constraint between the two phases of the traditional DAPFSP to form a just-in-time DAPFSP to minimize the maximum weighted tardiness cost. A variable neigh-

borhood search based on memory algorithm is proposed. From the above reviewed studies, it can be concluded that these scholars and researchers have extended and improved this problem, and most of them have an objective function based on minimizing the makespan, which provides a theoretical basis for the objective function in this paper.

Most scheduling problems are NP-hard and exact methods simply cannot find the optimal solution in reasonable computational time. Therefore, heuristic and meta-heuristic algorithms have been used to find the optimal solution and near-optimal solutions in reasonable computation time. Zhang et al. [11] pointed out that the DAPFSP is a typical NP-hard combinatorial optimization problem, and proposed an innovative three-dimensional matrix cube-based distribution estimation algorithm to address the difficulties of the proposed DAPFSP-T model. Likewise, Zhang et al. [12] also assumed this idea, and their team integrated the problem with two machine environments, distributed production and flexible assembly, which can assemble job processing into customized products, and proposed a mixed integer linear programming model to characterize the problem nature and solve the small-scale problem. An efficient modulo algorithm is further proposed. They affirm that the modal algorithm is an efficient algorithm to solve the DAPFSP, while some scholars prefer to use the greedy algorithm to solve the problem. For example, Pan et al. [13] solved a novel DAPFSP by proposing a mixed integer linear model, three construction heuristics, two variable neighborhood search methods, and an iterative greedy algorithm to obtain important problem-specific knowledge to improve the effectiveness of the algorithm. Similarly, Ochi et al. [14] studied the DAPFSP with the objective of minimizing completion time and, proposed an iterative greedy based approach called bounded search iterative greedy algorithm. Similarly, in order to coordinate production and transportation scheduling, Yang [15] proposed a novel DAPFSP-FABD model. The objective is to minimize the total cost of delivery and delay, and four heuristic algorithms, a variable neighborhood descent algorithm, and two iterative greedy algorithms are proposed. Liu et al. [16] proposed a memory algorithm for DAPFSP based on variable neighborhood search with the objective function of, i.e., makespan. And, the initialization based on NEH heuristic is applied to the product ordering. The neighborhood structure is introduced into VNS and used to perturb the job assignment of the factory and adjust the job order of the factory. Huang et al. [17] considered the DAPFSP with a total flow time criterion, and proposed an improved iterative greedy algorithm based on groupthink to solve the problem. It can be seen that the exact algorithm is no longer able to solve the problem, and this has promoted the development of heuristic and meta-heuristic algorithms, which can very definitely find the optimal solution as a solution to this scheduling problem.

More and more scholars are focusing more on the improvement of algorithmic aspects. For the no-wait flow shop scheduling problem that depends on time series, Hu et al. [18] proposed an enhanced differential evolutionary algorithm. Subsequently, Seidgar [19] took minimizing the weighted sum of expected completion time and average completion time as the solution objective and proposed four metaheuristic algorithms; namely, genetic algorithm, imperialistic competition algorithm, cloud theory-based simulated annealing, and adaptive differential evolution algorithm. Moreover, Li et al. [20] argued that the imperialist competition algorithm can solve the fuzzy distributed assembly flow shop scheduling problem. Furthermore, Al-Behadili et al. [21] proposed a multi-objective optimization model and particle swarm optimization solution method for the robust dynamic scheduling of permutation flow shops with uncertainty. Likewise, Zhang [22] emphasized that DAPFSP is a new generalization of the distributed displacement flow shop scheduling problem and the assembly flow shop scheduling problem, proposed an enhanced population-based metaheuristic genetic algorithm, and designed an effective crossover strategy based on local search to accelerate convergence. Li et al. [23] studied the minimization problem of DAPFSP and proposed a genetic algorithm with an enhanced crossover strategy and three different local searches. It is no coincidence that Mao et al. [24] advocated an improved discrete artificial bee colony algorithm to solve the DAPFSP with preventive maintenance operator, optimized by the completion time criterion. Moreover, a local search method

with insertion and exchange operators is used to generate adjacent solutions in the hiring bee phase and the bystander bee phase. Song and Lin [25] jointly proposed a genetic programming hyperheuristic algorithm to solve the DAPFSP with sequence-dependent setup times, minimizing the completion time as the objective. For the improvement of the algorithm, these scholars further study the genetic algorithm and particle swarm algorithm. This also provides a reference basis for our subsequent simulation experiments.

#### **3. Problem Description**

The traditional distributed assembly permutation flowshop scheduling problem (DAPFSP) can be divided into two phases; namely, the production phase and the assembly phase. In this section, the two-stage DAPFSP problem is extended to a three-stage DAPFSP, named the production stage, the transportation stage, and the assembly stage, and the problem is described as follows.

The job *j* ∈ {1, 2, ... , *n*} consists of operations *Oij*, *OTj* and *OAj*. In the production stage, there are *n* work processes and *f* identical plants, and job *j* ∈ {1, 2, ... , *n*} needs to be assigned to any plant *l* ∈ {1, 2, ... , *f* } for processing, and each plant is equipped with the same m machines, *M* = {*M*1, *M*2,..., *Mm*}, which correspond to the same assembly line shop for part processing. Operations *Oij* ∈ {*O*1*j*,*O*2*j*, ... ,*Omj*} are processed on machine *Mi*, *i* = 1, ... , *m* and require *pij* time units. Machine *Mij* can process at most one job at a time. The production and transport operations *OTj* will be executed on machine *MT* and take *pTj* time units. The rule is that in each assembly permutation flowshop, all jobs need to be processed on the same path in the order of machine *i*, which is equivalent to first on machine *M*1, then on *M*2. This goes on, until the end on machine *Mm*. Parts belonging to the same product must be processed continuously, and no partial parts are allowed to pass through. All jobs start timing at time *t* = 0 and on machine *i* = 1.

The subsequent phase is the transport stage, which means that after the completion of the first one in the first stage of the product operation, the transport machine collects all parts of the product and moves to the assembly stage.

In the assembly stage, there are *s* products and an assembly machine *MA*. Each product *r* ∈ {1, 2, ... ,*s*} has a defined part of the operation. The assembly operation *OAj* will be executed on the machine *MA* and uses *pAj* time units. The rule is that the assembly of a product can only start when the work contained in the product is completed and the assembly machine is idle. Continuous processing stages are part of the same product operation. The next product operation can only be performed after all jobs belonging to a product have been processed.

From this, it is clear that the objective function of the problem is to determine the allocation of products to plants, and the sequence of parts and products to each plant in order to minimize the completion time or maximum completion time for all products. Based on the above description of the problem, the minimization of the completion time is denoted as *Fm*|*nwt*|*C*max. To satisfy the above constraint, a feasible schedule *π* = {*π*1, *π*2,..., *πn*} is found for *n* jobs, such that the completion time *C*max(*π*) is minimized. This is the same where *C*max(*π*) is also equivalent to the time to complete the last job on the last machine. The Equation (1) is as follows.

$$\mathcal{C}\_{\text{max}}(\pi) = \sum\_{k=2}^{n} D\_{\pi\_{k-1}, \pi\_k} + \sum\_{j=1}^{m} P\_{\pi\_{n,j}} \tag{1}$$

In the transportation phase, it is necessary to ensure that the completion time distance between two adjacent jobs is determined by the processing time of the two processes, independent of the other jobs in the arrangement. Therefore, the completion time distance

is defined between each pair of jobs. The completion time distance between two adjacent job processes *πk*−<sup>1</sup> and *π<sup>k</sup>* is calculated as

$$\mathbb{C}\_{\pi\_{k-1},\pi\_k} = \max\_{1 \le i \le m} \left\{ \sum\_{j=1}^i P\_{\pi\_{k-1},j} - \sum\_{j=2}^i P\_{\pi\_k,j-1} \right\}, k = 2, 3, \dots, n \tag{2}$$

In order to build a mathematical model based on the above description, the constraints, parameters and variables are as follows.

Subject to,

$$\sum\_{j=0, j\neq k}^{n} X\_{j,k} = 1, k \in \{1, 2, \dots, n\} \tag{3}$$

$$\sum\_{\substack{j,k=0, k\neq j}}^n X\_{j,k} = 1, j \in \{1, 2, \dots, n\} \tag{4}$$

$$\sum\_{k=1}^{n} X\_{0,k} = f \tag{5}$$

$$\sum\_{j=1}^{n} X\_{j,0} = f \tag{6}$$

$$X\_{j,k} + X\_{k,j} \lessapprox 1, k \in \{1, 2, \dots, n-1\}, k > j \tag{7}$$

$$\sum\_{j=0}^{z\_{r-1}} \sum\_{k=z\_{r-1}+1}^{z\_r} X\_{j,k} + \sum\_{j=z\_r+1}^{z\_r} \sum\_{k=z\_{r-1}+1}^{z\_r} X\_{j,k} = 1 \tag{8}$$

$$r \in \{1, 2, \dots, \cdot, s\}, z\_r = \sum\_{t=0}^r n\_t \tag{9}$$

$$\mathcal{C}\_{i,j} \geqslant \mathcal{C}\_{i-1,j} + p\_{i,j}, i \in \{1, 2, \dots, m\}, j \in \{1, 2, \dots, n\} \tag{10}$$

$$\begin{array}{l} \mathbf{C}\_{i,j} \rhd \mathbf{C}\_{i-1,j} + p\_{i,j} + \binom{\mathbf{X}\_{j,k} - 1}{\mathbf{S}} \mathbf{g} \\ \dot{\mathbf{n}} \in \{1, 2, \dots, m\}, j \in \{1, 2, \dots, n\}, k \in \{0, 1, 2, \dots, n\}, j \neq k \end{array} \tag{11}$$

$$\sum\_{t=1, t \neq r}^{s} \mathbf{y}\_{r,t} = 1, r \in \{1, 2, \dots, \cdot, s\} \tag{12}$$

$$\mathcal{Y}\_{\mathbf{r},t} + \mathcal{Y}\_{\mathbf{t},\mathbf{r}} \lesssim 1,\\
\mathbf{r} \in \{1, 2, \dots, s - 1\},\\
\mathbf{t} > \mathbf{r} \tag{13}$$

$$\begin{array}{l} \mathbb{C}\_{r} \geqslant \mathbb{C}\_{m,j} + p\_{r} \\ j \in \{z\_{r-1} + 1, \dots, z\_{r}\}, r \in \{1, 2, \dots, \dots, s\} \end{array} \tag{14}$$

$$\begin{array}{l} \mathbb{C}\_{r} \geqslant \mathbb{C}\_{t} + p\_{r} + (\mathbb{Y}\_{r,t} - 1)\mathbb{g} \\ r \in \{1, \cdots, s\}, t \in \{0, \cdots, s\}, r \neq t \end{array} \tag{15}$$

$$TF \geqslant \sum\_{r=1}^{s} \mathbb{C}\_r \tag{16}$$

$$X\_{j,k} \in \{0, 1\}, j \in \{0, \dots, n\}, k \in \{0, \dots, n\}, j \neq k \tag{17}$$

$$\forall \, \mathbf{y}\_{r,t} \in \{0, 1\}, r \in \{0, \dots, \cdot, \text{s}\}, t \in \{0, \dots, \cdot, \text{s}\}, r \neq t \tag{18}$$

$$C\_{i,j} \gg 0, i \in \{1, 2, \dots, m\}, j \in \{1, 2, \dots, n\} \tag{19}$$

$$\mathbb{C}\_r \ni 0, r \in \{1, 2, \cdots, s\} \tag{20}$$

where *Xjk* represents a binary variable. It is equal to 1 if job *j* is the predecessor of job *k*. Otherwise, *Xjk* = 0. The same is true for *Yrt*. The product *r* is equal to 1 if it is the predecessor of the product *t*. Otherwise, *Yrt* = 0. *Cij* is the completion time of job *j* on machine *i*, and *Cr* is the completion time of product *r*. *pr* is the assembly time of product *r*, and *zr* is the total number of jobs contained in the first *r* products in the product sequence. In addition, *g* is a large enough positive number. *rt* is the number of jobs contained in product *t*. It is worth noting that *C*0,*<sup>j</sup>* = *Ci*,0 = *C*<sup>0</sup> = 0.

Equation (1) shows that the goal of the mathematical model is to minimize the completion time. The constraint sets (3) and (4) indicate that each job must have only one preceding job and one following job. Constraint sets (5) and (6) show that the preceding and following jobs of virtual job 0 are forced to be executed *f* times. Constraint set (7) ensures that a job cannot be both a predecessor and a successor of another job. Constraint set (8) ensures that jobs belonging to the same product are processed consecutively. Constraint set (9) indicates that job *j* can be processed on machine *i* − 1 only after it is processed on machine *i*. Constraint set (10) indicates that job *j* can only be processed on machine *i* after completing the processing of job *k* if job *j* is a successor to job *k*, where *g* is positive and has a sufficiently large scale volume. The constraint sets (11) and (12) ensure that each product must have only one preorder and one subsequent job. Constraint set (13) ensures that a product cannot be both a preorder and a successor of another product. Constraint set (14) ensures that each product cannot enter the assembly stage before the production of the part of the job it contains is complete. Constraint set (15) shows that if product *t* is a preceding job of product *r*, then product *r* cannot start the assembly job until product *t* is completed. Constraint set (16) defines the minimum completion time. Constraint sets (17) to (20) define the domains of the decision variables.

The model uses a sequence-based variable and a set of (min{ *f* ,*s*} + 1) virtual parts. The sequence starts with a virtual part and ends with a virtual part. These virtual parts divide all parts into min{ *f* ,*s*} subsequences, each corresponding to a plant. Parts belonging to the same product are never separated. Thus, the product sequence is implicitly included in the part sequence. For example, *s* = 4, *n* = 8, *m* = 2, *f* = 2. Table 1 shows the machining times of parts and products.


**Table 1.** Machining time for parts and products.

One of the feasible solutions is *x*<sup>04</sup> = *x*<sup>43</sup> = *x*<sup>35</sup> = *x*<sup>50</sup> = *x*<sup>01</sup> = *x*<sup>12</sup> = *x*<sup>28</sup> = *x*<sup>86</sup> = *x*<sup>67</sup> = *x*<sup>70</sup> = 1. The sequence of parts is {0, 4, 3, 5, 0, 1, 2, 8, 6, 7, 0}, consisting of the subsequences {4, 3, 5} and {1, 2, 8, 6, 7}.

#### **4. Materials and Methods**

To improve the local search capability of the algorithm, this section introduces the mathematical representation of feasible solutions, classification of populations, selection mechanism, domain structure and Gaussian learning strategy and elite learning strategy. A hybrid bat algorithm adapted to the three-stage DAPFSP problem is proposed.

#### *4.1. Mathematical Expression of the Feasible Solution of the Three-Stage DAPFSP*

The representation of feasible solutions is a crucial and important step in every metaheuristic approach. The addition of a transportation process to the DAPFSP solved in this paper. The result is that the original expression for the process permutation of the PFSP is not adequate for this scheduling problem. Therefore, in this section, the feasible solution is divided into a sequence of products and a sequence of *s* processes based on the model of the three-stage DAPFSP model and the properties of the bat algorithm. Each sequence of these s processes corresponds to a product. The processes of a part can be viewed as the sequence of product assembly on the assembly flowshop.

In the production and transportation stages, parts need to be allocated one by one according to the product order. When allocating parts, their processes are sequentially arranged to the corresponding factories according to the process order of the product, and a minimum manufacturing cycle is obtained. Only after all the processes included in the current part are completed can the next part be produced. When a product has completed all processes, the product can proceed to the assembly stage for the installation and assembly stages. As an example, suppose a set of sequences *f* is used to represent a feasible solution, and each plant has one and only one solution. Each sequence is an ordering of all the parts assigned to that plant, indicating the order in which the parts enter the flow shop. Thus, the order of processed products in the assembly stage is also implied, and the feasible solution can be expressed as *π* = *π*1, *π*2,..., *π<sup>f</sup>* , where

*π<sup>k</sup>* = *πk*,1, *πk*,2,..., *πk*,*η<sup>k</sup>* , *k* = 1, 2, . . . , *f* is the sequence associated with sequence associated with plant *k* and *η<sup>k</sup>* is the total number of parts assigned to plant *k*. For the example in Section 3, the representation corresponding to this feasible solution is *π* = {*π*1, *π*2} and *π*<sup>1</sup> = {4, 3, 5}, *π*<sup>2</sup> = {1, 2, 8, 6, 7}.

In the feasible solution of the original DAPFSP problem, the first stage represents the arrangement of all *N* processes, and the second stage represents the plants to which these *N* processes are assigned respectively. In this section, the concept of product priority is introduced to determine the order of product assembly, and a novel feasible solution of the three-stage DAPFSP is proposed. Specifically, the feasible solution in the first stage is expressed as the processing priority of the product, and the parts assembly of each product is equal to the number of processes that constitute that product. The specific steps are, firstly, calculating the total processing time of the product at each stage and randomly initializing the sequence of products. Second, the first two products in the initial product sequence are selected for evaluation, and whichever sequence is better is chosen as the current sequence. Again, this is done by placing it among the product sequences that are already scheduled properly. Finally, an optimal schedule can be obtained to select the current best sequence for the next iteration of the process. The second stage is represented as the processing priority of all processes, where the processes belonging to the same product are arranged in a random order. The smaller the number of its columns, the higher the priority of the element values; the third stage is the plant priority assignment vector, where each part is represented as the plant assigned to the corresponding process. The specific steps are to assign a part of the processes to the factories one by one in order, then select a process and find the minimum completion time by placing it after the last process in each factory. The best part assignment is selected for the next iteration.

Table 2 describes how the three-stage DAPFSP representation is decoded into a schedule. Suppose there are 3 products (*s* = 3) consisting of 9 jobs (*n* = 9) processed in 3 plants (*f* = 3). The mathematical representation is *P*<sup>1</sup> = {3, 4, 6}, *P*<sup>2</sup> = {1, 2, 8, 9}, *P*<sup>3</sup> = {5, 7}. According to the last two stages, called job sequence and factory assignment, it can be observed that jobs 5, 4, and 8 are assigned to the first factory and are processed in the order of 5, 4, and 8 according to the priority specified in layer 2. Then, it is instantly obtained that the partial scheduling of the factories can be decoded as *π*<sup>1</sup> = {5, 4, 8}. Similarly, the partial scheduling of the remaining two plants can be decoded as *π*<sup>2</sup> = {2, 7, 9} and *π*<sup>3</sup> = {6, 3, 1}, respectively. According to the values of the first stage specifying the product assembly

priority, the value of 1 in the table corresponds to jobs 5 and 7 belonging to product 3. It can be concluded that after the processing stage, product 3 should be assembled first. In a similar way, it can be concluded that product 1 is assembled for the second time immediately after product 2.

**Table 2.** Coding scheme.


In the context of the traditional representation of feasible solutions, the order of assembly of all products follows a first-come, first-served rule. This results in products being ranked in the processing stage in ascending order of the total completion time of their component operations. In contrast, the priority of products in the assembly phase is determined by the order of arrival at the assembly, and there is no way to specify it. This observation is confirmed by a study of the two-stage assembly scheduling problem by Tozkapan [26]. The efficiency may be higher if the sequence of operations in the processing phase is different from the sequence of operations in the assembly phase under the total delay objective. With a very late delivery date, for example, the assembly of the product may be delayed, even if the entire processing phase of the product's component operations has been completed.

#### *4.2. Population Classification*

The optimization process is divided into two stages. The first stage is that when the individual bat is in a better search position and is close to the optimal solution, its loudness and pulse emission rate reach the best state. This process is called the search phase, its population is called the search-type population. The second stage is the population of bats in a disadvantaged search position; that is, the capture population, so two populations with different functions are obtained. After each iteration, by updating the search direction and step length, the loudness and pulse emission rate of the bat algorithm are improved to find the optimal solution. By adjusting the weights and deviations of the network, the difference between the average value and the standard deviation of the bat algorithm can be minimized.

#### 4.2.1. Search-Type Population

The combination of the back propagation (BP) algorithm based on mean square error (MSE) and gradient descent method minimizes the mean square error [27]. The formula based on the mean square error measurement is shown in (21):

$$\text{MSE}(d, y) = \frac{1}{N} \sum\_{n=1}^{N} (d(n) - y(n))^2 \tag{21}$$

where *d*(*n*) represents the *n*-th element of the required signal, and *y*(*n*) represents the nth actual output. In the training process, the weight vector in the (*t* + 1) iteration is updated by (23). Where *α* is the learning rate or step length, W*<sup>t</sup>* is the weight vector of the previous iteration and *gt* is the gradient vector, which can be calculated by (22) and (23).

$$w\_{t+1} = w\_t - \alpha \mathfrak{g}\_t \tag{22}$$

$$\mathbf{g}\_t = \left. \frac{\partial \varepsilon}{\partial w} \right|\_{w = w\_l} = \left[ \frac{\partial \varepsilon}{\partial w\_{11}} \dots \frac{\partial \varepsilon}{\partial w\_{ij}} \dots \frac{\partial \varepsilon}{\partial w\_{nn}} \right] \tag{2.3}$$

In (23), *e* is the MSE error output in the *t*-th step of the training process and *∂*e/*∂*w is derived from the MSE error on each element of the *w* vector. The weight is expressed as *wt*+<sup>1</sup> = *wt* − *αtgt*, *α<sup>t</sup>* will be adjusted to an appropriate value to achieve better convergence. At this time, suppose the probability pops of the search-type population, its loudness and pulse emission rate are updated to:

$$A\_i^{t+1} = \mathfrak{a}\_t \times A\_i^t \tag{24}$$

$$r\_i^{t+1} = r\_i^o \left(1 - e^{-wt}\right) \tag{25}$$

#### 4.2.2. Captive Population

As a branch of the gradient descent method, the conjugate gradient method (CG) is applied to nonlinear unconstrained optimization problems. This method has strong local and global convergence [26–28]. According to the cloud computing resource scheduling problem, the CG method uses (26) to generate a weight sequence:

$$w\_{t+1}' = w\_t' + a\_t'd\tag{26}$$

where *α <sup>t</sup>* is the result of the line search method, which can be an exact line search or an inaccurate line search, which is expressed as a step size here. In (26), *dt* is in the descending direction, which is expressed as the search direction here, and its formula is shown in (27):

$$d\_{t+1} = -\mathcal{g}\_{t+1} + \beta\_t d\_t \tag{27}$$

In (27), *β<sup>t</sup>* is the conjugate parameter, *gt*+<sup>1</sup> is the gradient of the objective function concerning the weight at step *t* + 1, and represents the direction of the last step, and the first step is *d*<sup>0</sup> = −*g*0.

In the iterative process, the loudness *Ai* and the pulse emission rate *ri* will also change. As the bat moves closer to its prey, its loudness will decrease, and the pulse emission rate will increase. At this time, the probability of the catching population is *poph* and *pops* + *poph* = 1. In summary, the loudness and pulse emission rate of the bat algorithm are updated to:

$$A\_i^{t+1} = \pi\_t' \times A\_i^t \tag{28}$$

$$r\_i^{t+1} = r\_i^o \left(1 - e^{-w't}\right) \tag{29}$$

#### *4.3. Selection Mechanism*

When it comes to the bat algorithm as a learning algorithm, each individual in the population needs to learn from the individually optimal and globally optimal individuals in order to find the optimal solution. In the field of multi-objective optimization, there is no single optimal solution, and the optimal selection mechanism needs to be redesigned. Which strategy is used to update the individual optimal *pbest* and the global optimal *gbest* respectively, thus connecting the decision variable space and the objective space, is particularly important for weighing the exploration and exploitation capabilities of the algorithm evolution process. In this paper, we consider fully mining individuals with better diversity and convergence to assist populations to jump out of the local Pareto frontier. Three guides are selected from the searching and capturing populations to guide the entire individual bat flight, and how the guides are selected from the two populations is specified below.

Unlike the traditional bat algorithm, this algorithm first generates a set of *θ* with *ω* solutions and performs a local search process for the optimal solution in the set *θ*. Subsequently, IHBA enters the iterative phase. In the iteration, the solution *ϕ* is selected using the selection mechanism in the set *θ*. The solution *ϕ* is applied to the search phase and two partial sequences can be obtained, one for the solution *ϕSearch* with the products and processes removed, and the other for the remaining part of the solution *ϕCaptive*. The local search for the products is applied to *ϕCaptive*, and the suboptimal solution *ϕ Captive* is

generated. IHBA again reinserts the removed products and jobs into *ϕ Captive*. This process is called the capture phase, and a complete solution *ϕ* is regenerated. A local search is performed on *ϕ* to generate *ϕ* . Finally, acceptance criteria are used to determine whether the new solution *ϕ* updates the set.

Therefore, in scheduling problems, the selection mechanism is to select a solution from the solution set *θ* at the beginning of each iteration, and the selected solution is noted as *ϕ*, and proceeds to the subsequent stages such as search and capture. In this scheduling problem, the minimum completion time and the number of iterations are often used as selection factors. So, the proposed selection mechanism consists of two random options, both of which have a 50% probability of being selected. That is, two solutions are randomly selected in the set *θ* and the one with the lower minimum completion time is chosen as *ϕ*. Alternatively, two solutions are randomly selected in *θ*, and the one with the lower number of iterations is chosen.

The main focus of *pbest*, *gbest* and osd wizard selection mechanisms is on the target space, implying how to fully exploit the candidate solutions with good convergence and diversity in the target space. These wizards serve as a bridge between the target space and the decision variable space, and the new wizard selection mechanism can weigh the ability of exploration and exploitation of the decision variable space in an integrated way.

#### *4.4. Three Neighborhood Structures Based on Insert Operator*

To solve the three-stage DAPFSP problem, this paper introduces three neighborhood structures based on the INSERT operator, which constitutes the three-stage bat algorithm (3SBA). The insert operator is considered by many scholars as a superior performance operator, which is described as follows.


A process is randomly selected from the representation of the solution and assigned to another randomly selected plant. For this plant, all possible positions where jobs can be inserted are considered and the part of the job sequence with the earliest completion time is selected. It is worth noting that the sequence of jobs assigned to other plants remains unchanged during this process.

• Job Sequence-based Neighborhood (JSBN)

A plant is randomly selected and its sequence of jobs is extracted from the solution representation. Next, a job is randomly selected from this partial sequence, and another location is randomly inserted. Finally, the jobs in this partial sequence along with their product priorities and plants are placed in an orderly manner in the distribution solution representation, which is partially the location belonging to that plant.

#### *4.5. Local Search Method Based on Variable Neighborhood Structure*

Based on the above three neighborhoods named PBN, FBN and JSBN, this section proposes a local search method with variable neighborhood search. Equation *n*(*s*) denotes that the neighborhood structure and other parameters are consistent with the original bat algorithm, implying that the maximum number of iterations for each neighborhood at the same frequency is noted as *Tmax* and the frequency *f* is initialized as *f*0. In 3SBA, all three neighborhoods are iterated sequentially, avoiding premature convergence throughout the search process.

Variable neighborhood descent (VND) is the simplest variant of variable neighborhood search [29]. It has been shown to be effective in solving the flow shop [30] and distributed scheduling problems [31]. VND searches for different neighborhood structures from minimum to maximum. First, starting from the initial solution, VND searches the first neighborhood to find the local optimum. Then, while searching the second region, if a better local optimum solution is found, the algorithm returns to the first neighborhood. Otherwise, it continues searching the third neighborhood and so on. The algorithm does not end until a local optimum is found with respect to all neighborhoods. Therefore, to facilitate the differentiation of the VNDs in this paper, two versions were developed based on the characteristics and features of the bat algorithm. The first one uses three neighborhood structures, named VND3BA, and the second one uses two neighborhood structures, named VND2BA, respectively. A VND that best fits the bat algorithm and the scheduling problem is selected by calibration through simulation experiments.

#### 4.5.1. VND3BA Local Search Method Based on Three Neighborhood Structures

The first method is applied to each product contained in the plant. After removing each part of the product in succession, it is tested at all possible locations within the product. If a better manufacturing span is found, the part sequence is updated. This method runs until all parts contained in the plant have been taken into account and no further improvements are found; then, it can be stopped and ended. It can be named as Local Search Inside Product() and Part Local Search Inside Factory(). The former mainly searches individual products, while the latter searches the whole factory by calling the former, each in its own way.

The second method is a local search of the products within the factory. The principle of this method is to remove the products one by one from the beginning to the end and try to find the best position for them, equivalent to the optimal solution of this problem. If, in the process, a better manufacturing span is found, the sequence of parts in the plant is updated. The process is repeated until no products are found that can be improved.

The third method focuses on finding the most suitable location for a product that crosses plants. First, all plants are checked and the plant with the maximum completion time *p* is found. Then, the first product is removed and is tested in all possible locations in the remaining plants. If the generated best completion time is lower than the original completion time *Cp*, this product is inserted into the current location. The principle is that the algorithm iterates by finding the plant with the maximum completion time again, until a local optimum is found.

#### 4.5.2. VND2BA Local Search Method Based on Two Neighborhood Structures

Since VND3BA considers the neighborhood of parts and the neighborhood of products separately, some available feasible solutions may be missed. Therefore, in this section, two local search algorithms are proposed to consider these two neighborhoods in an integrated manner.

The first approach explores the interior of the plant by extracting the range of products and finding the best location of the products within the plant. This step is to extend the scope of the local search. This method also adaptively finds the best part order for the product when inserting it.

The second approach exploits the local search capability of the bat algorithm itself to explore the optimal solution across plants. Therefore, it can be viewed as a combination of two processes named ProductLocalSearchBetweenFactories() and LocalSearchInsideProduct().

#### *4.6. Gaussian Learning Strategy and Elite Learning Strategy*

Gaussian learning strategy (GLS) [32] and elite learning strategy (ELS) [33] can assist individuals of bats to jump out of the local frontier and improve the local search ability of the algorithm. If the feasible solution is not updated many times, the whole bat population is most likely to fall into local optimum, and then only the GLS reset strategy can be executed, which is given by

$$\left|\mathbf{x}\_{i}^{t+1} \sim N\left(\frac{\mathbf{x}\_{gbest} - \mathbf{x}\_{gbest,i}}{2}, \left|\mathbf{x}\_{gbest} - \mathbf{x}\_{gbest,i}\right|\right)\right.\tag{30}$$

The collaborative learning mechanism of *pbest* and *bbest* ensures that the exploration and exploitation capability of the entire bat population is enhanced, allowing the algorithm to quickly jump out of the local frontier to approximate the true frontier.

Although the use of GLS can assist individual bats to jump out of the local frontier, this strategy alone still cannot meet the requirements of the bat algorithm for solving the DAPFSP. To further improve the performance of the bat algorithm, ELS is introduced to solve the scheduling problem with the following mathematical expressions.

$$E\_i(j) \sim E\_i(j) + (\mathbf{x}\_{ub}(j) - \mathbf{x}\_{lb}(j))N(0, 1) \tag{31}$$

where *xub*(*j*) and *xlb*(*j*) represent the upper and lower bounds of the *j*-th dimensional decision variable, respectively, and *N*(0, 1) represents the random number with a mean value of 0 and a variance of 1.

#### **5. Simulation Results**

#### *5.1. Test Questions and Parameter Settings*

In this section, the proposed IHBA is compared with the competitive memetic algorithm (CMA) [34], biogeography-based optimization algorithm (BBO) [35], estimation of distribution algorithm ( EDA) [36], genetic algorithm [19,21], and particle swarm optimization algorithm [20], and only one of these algorithms was selected for comparison from the literature [37]. So far, the three DIWO algorithms have proved to be the state-of-the-art algorithms for the considered problem with the parameters shown in Table 3. All simulation experiments were implemented on the MATLAB r2020a platform. All programs were executed on an AMD A8-7100 Radeon R5@1.80 GHz and 12.0 GB RAM in Windows 10 operating system.


**Table 3.** Algorithms and Parameters.


**Table 3.** *Cont.*

To calibrate the proposed IHBA, a total of 810 instances were randomly generated according to the literature [5,34] written by Hatami et al. These instances can be found in http://soa.iti.es/problem-instances, accessed on 1 September 2021. These instances combine four factors, including the number of jobs (*n*), the number of machines (*m*), the number of products (*s*), and the number of plants (*f*). Each factor has three levels, *n* ∈ {100, 200, 500}, *m* ∈ {5, 10, 20}, *f* ∈ {4, 6, 8} and *s* ∈ {30, 40, 50}. Thus, a total of 3<sup>4</sup> = 81 factor combinations are studied. The processing time in the production stage is randomly generated in the range [1, 99] with uniform distribution. The operation time of the transportation stage is obtained in the range [l, 99l] with uniform distribution. The assembly operation time for the assembly stage is obtained in the range between *nl* and 99*nl* with uniform distribution.

#### *5.2. The Results and Discussion*

Relative percentage deviations are the absolute deviations of a given measurement as a percentage of the mean, and can only be used to measure the deviation of a single measurement from the mean. In order to compare the efficiency between algorithms, the calculation of relative percentage deviations (RPD) is considered to compare and analyze these metaheuristic algorithms. Its formula is shown in (32).

$$RPD = \frac{\mathbb{C}\_{\text{max}}(\pi) - \mathbb{C}\_{\text{max}}(\pi)\_{\text{best}}}{\mathbb{C}\_{\text{max}}(\pi)\_{\text{best}}} \times 100\% \tag{32}$$

where *Cmax*(*π*) is the total completion time of the current algorithm, and *Cmax*(*π*)*best* is the minimum value of *Cmax*(*π*) among all algorithms to be compared.

In this paper, the termination condition of the iteration is set to a maximum CPU runtime equal to *C* × *m* × *n* milliseconds, and the results are calculated for *C* = 20, 40, and 60. Then, the results obtained from the calculation are compared with other algorithms. The methods used are the analysis of variance (ANOVA), which is widely used in parametric statistical procedures, and the Friedman test and Wilcoxon paired signed test, which are commonly used in nonparametric statistical procedures. These methods have been widely used and promoted in the recent scheduling literature.

#### 5.2.1. Comparative Analysis at *C* = 20

When *C* = 20, the values of RPD for various heuristic algorithms are shown in Table 4. It can be seen from the bold figures that the value of PRD of the proposed algorithm IHBA in this paper is smaller than that of other algorithms in the vast majority of combinatorial solutions, including the recognized TDIWO algorithm and the original bat algorithm. Only in these 5 sets of data (1) *m* = 5, *n* = 100,*s* = 30 and *f* = 4, (2) *m* = 5, *n* = 200,*s* = 40

and *f* = 4, (3) *m* = 10, *n* = 100,*s* = 50 and *f* = 8, (4) *m* = 15, *n* = 200,*s* = 50 and *f* = 4, (5) *m* = 15, *n* = 200,*s* = 40 and *f* = 8, IHBA is 5–10% smaller than the other algorithms. Meanwhile, IHBA not only obtains the lowest relative percentage deviation value, but also the smallest average RPD value, which is 21.15%, with an error of about 0.83% from the optimal value. This is good proof that the IHBA algorithm has a strong advantage and superiority. It can also be seen that the GA and PSO performance is also very good, with average RPD values of 32.09% and 31.52%. Although TDIWO is not satisfactory in this comparison, it is significantly better than BBO and EDA. The most surprising is the original bat algorithm, with a result of 27.69, which is second only to the algorithm proposed in this paper. This also shows that the bat algorithm itself has better performance and has stronger search capability.

**Table 4.** The Average RPD Values at *C* = 20. It can be seen from the bold figures that the value of PRD of the proposed algorithm IHBA in this paper is smaller than that of other algorithms in the vast majority of combinatorial solutions, including the recognized TDIWO algorithm and the original bat algorithm.


**Table 4.** *Cont.*


Then, the results of the descriptive statistics of the Friedman test were calculated for the case of *C* = 20, as shown in Table 5, where *N* is the number of test cases and takes the value of 4050. The minimum value of each item derived in the algorithm is marked in bold. The *p*-value calculated for the test statistic is 0 when the significance level *α* = 0.05, which is less than or equal to 0.05. This indicates that there are significant differences between the algorithms used for comparison. The following conclusions can be clearly seen through Table 5. The rank of IHBA is only 2.60, which is close to one-half of CMA, BBO and TDIWO, while it is smaller than the ranks of GA, PSO and the original bat algorithm by 1.39, 1.09 and 0.78, respectively. Although IHBA does not take the best value in terms of maximum value, it means that standard deviation and minimum value are the smallest among all algorithms. The comparison between the groups also shows that the ranking of EDA is even more than three times that of IHBA, indicating the worst performance of the algorithm, which echoes the results and analysis in Table 4. In addition, it can also be seen by the mean and standard deviation that the three algorithms of DIWO are very close to the proposed algorithm.

**Table 5.** The Descriptive Statistics Achieved by the Friedman Test at *C*=20 and *α*=0.05. The bold figures indicate the comparison between algorithms.


Next, the non-parametric test for paired samples is the Wilcoxon paired signed test method used. This method was developed based on the signed test for paired observations, and is more effective than the traditional test with positive and negative signs alone. The observed data are generally considered to be significantly different when *p* is less than 5%. When *p* is greater than or equal to 5%, the difference in the data is considered insignificant. The results of the Wilcoxon paired signed test for this paper are given in Table 6, assuming that there is no difference between each pair of algorithms (*α* = 0.05). *R*+ in the table represents positive differences and *R*− represents negative differences. In the results for each pair of algorithms, the sum of *R*+, *R*− and bonding values is equal to *n* in Table 5. It can also be seen in Table 6 that the *p*-values in the last column are equal to 0 and all are less than 0.001, indicating that the differences between the algorithms compared are statistically significant with respect to 0. In other words, there is a significant difference between IHBA and the other algorithms.

**Table 6.** The Wilcoxon paired signed test result at *C* = 20 and *α* = 0.05.


Finally, in order to make the observations more visual and relevant, this section analyzes and interprets the results of the ANOVA analysis for the eight algorithms. The methods used are the mean plot and the minimum difference method interval, and Figure 1 shows the plotted plots with 95% confidence level. It can be very clearly seen that the best performance is IHBA, which beats the other comparison algorithms by a more significant margin. The bat algorithm also has some advantages, while GA and PSO continue to perform consistently, and the EDA algorithm has the worst performance, followed by BBO.

**Figure 1.** Mean plot and 95% LSD interval of the bat algorithm at *C* = 20.

#### 5.2.2. Comparative Analysis at *C* = 40

The lower the RPD, the better the performance of the algorithms. When C takes the value of 40, the average RPD of these eight algorithms is shown in Table 7. As can be seen in the last row, the average value of IHBA is 22.73%, which is significantly lower than the average value of other algorithms including the original bat algorithm. It shows that the IHBA algorithm is more efficient than the other metaheuristic algorithms, while BA performs second, with an average RPD of 23.57%. Compared with Table 7, it can be concluded that the RPD value of IHBA changes from the original 21.15 to the current 22.73, as the value of C is taken to increase, which shows that the value of RPD of IHBA increases with the increase of the number of iterations and time. In addition, PSO with a mean value of 35.76% is significantly better than GA and the other four heuristics. Meanwhile, EDA has the largest value and its performance is the worst, which is also consistent with the results in Table 4.

Table 8 describes the results of the Friedman test calculations. As can be seen at a glance from the rankings, IHBA has a mean rank of 3.06, which is the smallest value among these algorithms. The next best performer is TDIWO, which also has a better performance, but its standard deviation is 0.89, which is larger than CMA, GA, and PSO. It can be seen from the mean and standard deviation that IHBA has values of 0.79 and 0.76, respectively, which are also the smallest values among all algorithms. The performance of EDA is still the worst, with its average rank, mean, standard deviation, and maximum value all being at a disadvantage. The performance of EDA is still the worst, with its ranking, mean, standard deviation and maximum value being at a disadvantage. In particular, the mean value is as high as 32 times that of IHBA. The *p*-value in the last row is 0, which is less than 0.05, indicating that the difference between the algorithms has statistical significance. This difference is not due to chance sampling, but rather, the difference between the two groups of algorithms is significant.


**Table 7.** The Average RPD Values at *C* = 40. It can be seen from the bold figures that the value of PRD of the proposed algorithm IHBA in this paper is smaller than that of other algorithms in the vast majority of combinatorial solutions, including the recognized TDIWO algorithm and the original bat algorithm.

**Table 7.** *Cont.*


**Table 8.** The Descriptive Statistics Achieved by the Friedman Test at *C* = 40 and *α* = 0.05. The bold figures indicate the comparison between algorithms.


The results obtained from the Wilcoxon paired signed test are shown in Table 9. The *p*-values in the last column are all equal to 0, which is less than 0.05. It shows that the differences between these algorithms are statistically significant, and that there are significant differences between IHBA and the other algorithms. This also reconfirms the conclusion that *C* = 20. More definitely, it explains that IHBA is an effective group intelligence algorithm with outstanding performance.


**Table 9.** The Wilcoxon paired signed test result at *C* = 40 and *α* = 0.05.

Finally, in this section, EDA is excluded from the mean plot. The reason is that the algorithm is too intrusive and its performance is too poor. As can be seen in Figure 2, IHBA is once again stronger than the other metaheuristic algorithms by a margin. Once again, the variability between the algorithms is proven to be at a significant level. Since the LSD test has the highest sensitivity, Figure 2 verifies from the side that IHBA is an algorithm that performs well and can reach Pareto optimality.

**Figure 2.** Mean plot and 95% LSD interval of the bat algorithm at *C* = 40.

5.2.3. Comparative Analysis at *C* = 60

Table 10 presents the results of the calculations at C = 60. Table 10 shows the validity of the IHBA, which yielded an overall mean RPD value of 31.55%. The results of the Friedman test in Table 11, the results of the Wilcoxon paired sign test in Table 12, and the mean plot in Figure 3 again show that the proposed IHBA performed best in the comparison.


**Table 10.** The Average RPD Values at *C* = 60. It can be seen from the bold figures that the value of PRD of the proposed algorithm IHBA in this paper is smaller than that of other algorithms in the vast majority of combinatorial solutions, including the recognized TDIWO algorithm and the original bat algorithm.

**Table 10.** *Cont.*


**Table 11.** The Descriptive Statistics Achieved by the Friedman Test at *C* = 60 and *α* = 0.05. The bold figures indicate the comparison between algorithms.


**Table 12.** The Wilcoxon paired signed test result at *C* = 60 and *α* = 0.05.


**Figure 3.** Mean plot and 95% LSD interval of the bat algorithm at *C* = 60.

#### *5.3. VND Method*

In this section, the proposed VND3BA and VND2BA are compared with VNDH12, VNDH22, and VNDH32, proposed by Hatami et al. These are 3 different versions of the VND method using H12, H22 and H32 to generate the initial solution. The VND method uses product local search to improve the product sequence and part local search to improve the part sequence for each product. Tables 13 and 14 show the computational results and CPU time for the 5 VND algorithms.

As the table shows, VND3BA and VND2BA perform significantly better than VNDH12, VNDH22 and VNDH12 in terms of solution quality and computational effort. The best overall average RPI value of 15.233% for all instances obtained by the existing VND algorithm is more than three times higher than the best overall average RPI value obtained by VND3BA and VND2BA. For all values of *n*, *m*, *f* and *s*, the proposed VND algorithm produces better results than the existing VNDH12, VNDH22, and VNDH32. The VND method shows stable advantages without considering the different *n*, *m*, *f* , and *s* involved. The proposed VND algorithm is also very fast, with an overall average CPU time of 0.009 and 0.057 s for VND3BA and VND2BA, respectively, compared with an overall average CPU time of more than 5 s for the existing algorithms.


**Table 13.** Average RPD value of VND method.


**Table 14.** CPU time of VND method.

#### **6. Conclusions**

To address the challenge that existing heuristic algorithms and meta-heuristic algorithms cannot fundamentally solve the three-stage distributed assembly permutation flow shop optimization problem, this paper proposes a hybrid bat algorithm optimization algorithm based on variable neighborhood structure and two learning strategies to minimize the completion time of this problem. The algorithm designs a search-based and capturebased two populations to solve the difficult trade-off between convergence and diversity of the bat algorithm in solving the scheduling optimization problem. Moreover, by fully mining the population information, a new selection mechanism and a new velocity and location update strategy are designed to solve the difficult trade-off between exploration and exploitation when the bat algorithm solves the scheduling optimization problem. The Gaussian learning strategy and elite learning strategy are used to assist the whole population to jump out of the local optimal frontier. Simulation results show that IHBA can solve the optimization problem of three-stage distributed assembly permutation flow shop scheduling well, and performs better than existing algorithms in the literature.

In future research, the algorithm is extended to more complex multi-objective optimization problems to further improve the efficiency of iterative search, and the proposed algorithm is applied to other scheduling problems.

**Author Contributions:** Conceptualization, Y.W.; methodology, Y.W.; software, Y.W.; validation, Y.W. and J.Z.; formal analysis, Y.W.; investigation, Y.W.; resources, Y.W. and J.Z.; data curation, Y.W.; writing—original draft preparation, Y.W.; visualization, Y.W. and J.Z.; supervision, Y.W. and J.Z.; project administration, J.Z.; funding acquisition, J.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Fundamental Research Funds for the Central Universities under Grant No. CUSF-DH-D-2018050 (18D310804).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The benchmark of instances for the permutation flowshop scheduling problem can be found here: http://soa.iti.es/problem-instances, accessed on 1 September 2021.

**Acknowledgments:** We gratefully acknowledge the anonymous reviewers for their insightful comments on the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **A Multi-Criteria Assessment of Manufacturing Cell Performance Using the AHP Method**

**Zuzana Soltysova \*, Vladimir Modrak and Julia Nazarejova**

Faculty of Manufacturing Technologies, Technical University of Kosice, Bayerova 1, 080 01 Presov, Slovakia; vladimir.modrak@tuke.sk (V.M.); julia.nazarejova@tuke.sk (J.N.)

**\*** Correspondence: zuzana.soltysova@tuke.sk

**Abstract:** Research of manufacturing cell design problems is still pertinent today, because new manufacturing strategies, such as mass customization, call for further improvement of the fundamental performance of cellular manufacturing systems. The main scope of this article is to find the optimal cell design(s) from alternative design(s) by multi-criteria evaluation. For this purpose, alternative design solutions are mutually compared by using the selected performance criteria, namely operational complexity, production line balancing rate, and makespan. Then, multi-criteria decision analysis based on the analytic hierarchy process method is used to show that two more-cell solutions better satisfy the determined criteria of manufacturing cell design performance than three less-cell solutions. The novelty of this research approach refers to the use of the modification of Saaty's scale for the comparison of alternatives in pairs based on the objective assessment of the designs. Its benefit lies in the exactly enumerated values of the selected criteria, according to which the points from the mentioned scale are assigned to the alternatives.

**Keywords:** multi-criteria assessment; cell manufacturing design; operational complexity; makespan; production line balancing rate

#### **1. Introduction and Problem Description**

The smart manufacturing concept, as an important part of the Industry 4.0 strategy, opens new opportunities for producers to implement new platform-based business models by embracing cutting-edge technologies. Cellular manufacturing (CM) systems belong among six fundamental manufacturing systems, for which Industry 4.0 has been conceived [1]. Their goal is to complete jobs as swiftly as possible, make a wide variety of similar products, and produce as little waste as possible. An important objective of CM systems is to be as easily reconfigurable as possible [2–4]. If the number of potential configuration design solutions is generated, the role of the user is to define constraints on design and performance to obtain solutions meeting their requirements. Subsequently, it is needful to explore design options and configurations to select an optimal solution.

Operation management research often reflects CM problems in order to make manufacturing operations more efficient and productive. The cellular manufacturing method brings scattered processes together with compact cells usually arranged in U-shapes and can significantly improve the operation of batch production [5]. Batch production with a wide range of product types presents a crucial problem for layout designers since the parts move in batches from one process to another, and ready parts must wait for the remaining parts to complete processing before they move to the next stage. Because operation times are distinct one from another, it causes unbalanced machine utilization, scheduling problems, and possible late deliveries. Unbalanced machine utilization is frequently considered one of the main criteria in CM design optimization since ineffective utilization of machines can result in unprofitable production [6,7]. Another approach to optimize CM design is based on using appropriate scheduling techniques that aim to determine the actual assignment

**Citation:** Soltysova, Z.; Modrak, V.; Nazarejova, J. A Multi-Criteria Assessment of Manufacturing Cell Performance Using the AHP Method. *Appl. Sci.* **2022**, *12*, 854. https:// doi.org/10.3390/app12020854

Academic Editor: Aliaksandr Shaula

Received: 15 December 2021 Accepted: 13 January 2022 Published: 14 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of products or jobs to available operations in order to complete manufacturing orders on time [8]. These techniques are mostly developed using the minimum makespan criterion for CM design optimization [9]. Moreover, makespan minimization positively influences work-in-process inventories [10,11]. Taking into consideration the fact that CM systems have to be adaptable to changing market demands, their structural and operational complexity adequately increases as a consequence [12,13]. In this context, the cost of operational complexity is becoming a topical problem, and therefore, the trade-offs between the level of complexity and its added value are often discussed (see, e.g., [14–17]).

It is also worth mentioning that CM is not always the appropriate approach to take for certain scenarios [18]. This problem has been studied, e.g., by Flynn and Jacobs [19], who indicated that, through a well-organized job shop, it is possible to achieve at least as good a performance as the cellular layout with respect to several criteria such as work-in-process inventory levels and average flow times. For this reason, it is assumed in this study that a given manufacturing process satisfies basic criteria for a cellular layout, which include a high ratio of setup to process time, stable demand, unidirectional workflow within a cell, and a considerable level of material movement times between process departments [20].

The problem that is treated in the present work can be described as follows. The time required to complete all jobs, namely makespan, is one of the most frequent performance indicators for the job shop problem [21]. Considering that one-piece flow cellular manufacturing is the ultimate in lean production [22], CM is primarily aimed at reducing times within the production system. On the other hand, competitiveness through cost reduction in the design and implementation of production systems is an important and permanent task for process designers. In this order, identifying and monitoring cost items causes significant difficulties since there are various costs such as part holding cost at a facility, machine procurement cost, machine maintenance overhead cost, machine repair cost, production loss cost due to machine breakdown, machine operation cost, setup cost, tool consumption cost, inter-cell travel cost, intra-cell travel cost, etc. [23]. Therefore, from a practical viewpoint, it is reasonable to indicate the main portion of costs by using indirect indicators, namely production line balancing rate and minimization of operational complexity.

Taking into account the above-mentioned observations and factors influencing the efficiency of CM systems, our interest in this paper is to propose a multi-criteria decisionmaking approach for the selection of an optimal manufacturing cell from several alternative options. The three selected criteria, i.e., maximization of production line balancing rate, minimization of operational complexity of alternative CM designs, and makespan minimization, is employed to assess the layout design alternatives.

For the purpose of multi-criteria decision analysis, the analytic hierarchy process (AHP) method is applied. The AHP method is one of the most exploited multi-criteria decision-making tools and one of the most trusted decision-making methodologies. On the other hand, its disadvantage lies in possible flaws in the verbal scale often used in AHP pairwise comparisons. In the case of investigation of the given manufacturing problem, a numerical scale is used for the assessment of the criteria. Then, the results will not be influenced by personal opinions in considering and representing facts.

As testing examples, alternative CM designs along with input data from the existing representative case study is used.

#### **2. Literature Background**

The literature background conducts a review regarding the above-mentioned optimization criteria or factors influencing decision making in cellular manufacturing design.

Line balancing is considered an effective tool to optimize layout design and reduce product cycle times. The objectives of line balancing techniques in flow shop scheduling problems are usually focused on the reduction and/or redesign of workstations in order to minimize production costs, work in process in order to reduce storage space and bound capital, and minimization of makespan and flow times [24,25]. In general, line balancing techniques are divided into deterministic types, where all input parameters are known and not changed over time, and probabilistic types, which deal with parameter uncertainties [26]. An innovative stochastic line balancing method was proposed by Kottas-Lau et al. [27]. Their algorithm was developed for the purpose of achieving optimized total production costs and to allow a good level of line balancing. Among the other important deterministic line balancing techniques can be mentioned the largest candidate rule [28], the Kilbridge and Wester Column method [29], and the ranked positional weight method proposed by Hegelson and Bernie [30]. The latter method was developed to minimize idle times and the number of workstations. The specific type of these optimization tasks consists of large-scale line balancing problems that deal with uncertainty in line balancing. Hazr and Dolgui [31] proposed two optimization models that belong to these problems, by which it is possible to generate exact solutions of cellular manufacturing designs. Another optimization technique based on the identification of equilibrium between line efficiency and equipment cost was proposed by Gurevsky et al. [32]. It is also useful to note that their method supports decisions at the early design stage of production lines. As it is impossible to involve all pertinent works on the topic, comprehensive literature studies focused on the comparative evaluation of line balancing techniques overcome this drawback. Such review studies can be found, e.g., in [33–35].

Operational and structural complexity is another relevant factor affecting the efficiency of CM systems. According to Hon [36], the main reason for the investigation of manufacturing system complexity is to comprehend and control the behavior of such systems in order to make them more productive and predictive. Fredendall and Gabriel [37] argue that by measuring system complexity, managers can better identify problems in manufacturing systems that hinder production flow. Therefore, the determination of quantitative metrics of manufacturing system complexity, either static or dynamic, is one of the crucial elements in this effort [38]. There are several different approaches to the utilization of complex concepts in this application domain, but it has not yet been possible to find closed-form equations able to describe the dynamic behavior of manufacturing systems [39]. For this reason, available operational complexity measures reflect only selected facets of such systems. Manufacturing system complexity is often divided into structural and operational types. The second type of complexity measure, which is of interest in this study, is based on measuring uncertainties involved in manufacturing systems. This type of complexity is further divided into time independent and time dependent [40]. Zhang [41] analyzed the relationship between cellular manufacturing system complexity and utility in order to show that increasing complexity can be beneficial for manufacturers until it reaches a critical value. Beyond this critical value, the situation becomes the opposite. The mentioned system complexity consequences on design and managerial practice were originally introduced by Tainter [42] and are widely respected in many research communities. Therefore, managers might mitigate the negative aspects of complexity while managing its positive aspects, as complexity indirectly influences the performance of manufacturing systems [43]. These arguments inspire us to employ operational complexity as one of the criteria in the decision-making procedure.

Makespan is commonly used as a criterion of performance measure in the design of cellular manufacturing systems, because the advantages of cellular manufacturing also include simplified planning and scheduling [44]. The scheduling problem in a cellular manufacturing system assumes that intercellular moves can be eliminated by duplicating machines, but it is usually very costly and therefore infeasible [45]. If duplicating machines is not a viable solution, then a volume limit can enhance the choice of the optimal routing of jobs. One method to make this choice is to minimize the makespan since this performance measure is the most frequent objective in flow shop problems [46].

#### **3. Methodological Framework**

This section aims to describe in a nutshell the set of methods and overall procedure in chronological order and help the reader understand the context under which the research was conducted. In its first steps, the criteria or factors influencing the efficiency of CM systems have been mostly identified based on empirical findings and general knowledge of CM design, namely makespan and production line balancing rate. Similarly, measurement methods to quantify operational complexity have been chosen. Subsequently, the AHP method divided into three steps was used for the evaluation of individual CM alternatives. Summarily, the methodological framework of this research consists of five steps, which is illustrated in the figure below (Figure 1).

**Figure 1.** Overall research procedure in chronological order (Sections 1–5).

*3.1. Indicators Used as Criteria to Assess Alternative CM Designs*

3.1.1. Production Line Balancing Rate

The production line balancing rate (PLB) represents a measure of the average length of time of every cycle time in the working procedure on the processing line. It is equivalent to minimizing the number of workstations with a certain takt time [47]. It is expressed as [48]:

$$PLB = \frac{\sum\_{j=1}^{n} t\_j}{m \ast \max(T\_i)} \text{ [\%]}\_r \tag{1}$$

where

*tj* stands for standard work time of the *j*-th job elements; *n* represents the number of the work elements; *m* represents the number of total lines in the production system; *Ti* represents the work time in the production line(s) (PL(s)); max(*Ti*) represents the biggest line operating time.

3.1.2. Operational Complexity Indicators

There are several potential operational complexity indicators that can be applied to measure the operational complexity properties of manufacturing systems from different or similar perspectives (see, e.g., [49–54]). Two of them that suit the available data characterizing the benchmarked alternatives are further introduced.

#### **Process Complexity Indicator**

This indicator, similar to the other concurrent complexity measures, is derived from Shannon's information theory [55]. The process complexity indicator (PCI) is specified for the quantification of manufacturing process complexity, taking into account the operational complexities of individual machines. The PCI indicator is enumerated by using equation [54]:

$$PCI = -\sum\_{i=1}^{M} \sum\_{j=1}^{P} \sum\_{k=1}^{O} p\_{ijk} \cdot \log\_2 p\_{ijk} \text{ [bits]},\tag{2}$$

where

*pijk* stands for the probability that part *j* is processed due to operation *k* by individual machine *i* according to scheduling order;

*O* is the number of operations according to parts production;

*P* is the number of parts produced in the manufacturing process;

*M* is the number of all machines of all types in the manufacturing process.

#### **Balanced Complexity Indicator**

The balanced complexity indicator (BCI) takes into account the rate of mutual differences between the individual complexities of machines. This indicator expresses the variability of the partial complexities of workstations/machines and calculates the deviation of partial machine complexities from their mean value. It is calculated using the following formula [54]:

$$BCI = \frac{\sum\_{i=1}^{N} MCI\_{i(\text{max})} - \sum\_{i=1}^{N} MCI\_{i(\text{min})}}{N} \text{ [bits]}.\tag{3}$$

where

*MCIi*(max) represents the first *N*-max complexity values; *MCIi*(min) represents the first *N*-min complexity values; *N* represents the number of max and min machine complexity values.

#### 3.1.3. Makespan Indicators

For the purpose of calculating makespan, the scheduling algorithm to minimize the completion of n-jobs of m-machines is used. As known, there are many different algorithms for the given purpose. In this work, the freely available online software is utilized [56], and the following input data for this algorithm are collected:

Processing times in minutes for each job, which are included in matrix *mxn* (Figure 2a); Number of transport batch (Figure 2b);

Transport batch sizes for each job (Figure 2c);

Sequences of individual jobs numbered by order (Figure 2d).

Its application in the first step requires that input data are presented in the form of a Microsoft Office Excel table and pasted into the input data window. A flowchart of how to generate makespan is shown in Figure 3.

As can be seen from Figure 3, makespans are enumerated assuming two scenarios. According to the first scenario, makespan is calculated for determined batch sizes, and for the second scenario, the one-piece flow (OPF) principle is applied, i.e., when the transport batch size for each job equals 1.

**Figure 2.** Table format for insertion of input data, (**a**) processing times in minutes for each job, which are included in matrix *mxn*, (**b**) number of transport batch, (**c**) transport batch sizes for each job, (**d**) sequences of individual jobs numbered by order.

**Figure 3.** Software flowcharts with acquired data.

#### *3.2. Description of AHP Method*

Prior to describing the AHP method with the modified Saaty scale for a comparison of design alternatives in pairs using a multi-criteria decision-making approach, the five selected related AHP approaches are reported in Table 1 in order to point out the differences in the research addressed in our paper.


**Table 1.** Comparison of existing studies based on usage of AHP method.

All of these AHP method applications are based on a subjective approach in mutual comparison of alternatives. On the contrary, the proposed procedure to find the most suitable cell design alternative(s) uses an objective approach of the pairwise comparison of cell design alternatives. In addition, it is necessary to select one of the possible methods in the application of the AHP method [63,64]. For the given purpose, the most accurate procedure seems to be that which consists of the following steps:

Step 1. Creation of hierarchical structure of AHP method. The hierarchical structure of the AHP method is created in the form of a diagram, where the criteria, sub-criteria, and CM alternatives are specified.

Step 2. Pairwise comparison of CM alternatives. All CM alternatives are pairwise benchmarked with respect to the criteria and sub-criteria. The pairwise comparison is provided in matrix form by comparing one CM alternative to another to determine the weights of importance. For this purpose, the proposed modified scale of relative importance is applied (see Section 4.3).

Step 3. Enumeration of priority vectors and aggregated results. The priorities are derived using the values of the principal right eigenvectors of the compared matrices. These priorities are expressed as absolute numbers bounded between 0 and 1, without units, and are calculated according to the so-called additive normalization method using the following simple procedure:

Sum each column values separately for each matrix, divide each element of the column with the sum of that column for each matrix, and compute the average of all elements in each row of all matrices to obtain the priority vector.

To obtain aggregated results, it is needed to summarize the determined priorities of all the individual indicators for each CM alternative. Then, the aggregated priorities are compared by ranking them in order from most to least important. Finally, the optimal CM alternative is selected.

#### **4. A Practical Example**

#### *4.1. Description of Manufacturing Cell Designs Alternatives*

This section aims to introduce five manufacturing cell alternatives for their mutual comparison in order to determine the optimal design. The CM alternative designs along with basic input data are taken from a case study by Yan and Irani [65]. The essential data from this case study include routing of parts (P) through machines (M), sequencing of parts by way of recording, operational times in minutes, and batch sizes for individual parts. The parts routings and operational times for all the parts are shown in Figure 4 (e.g., part P1 is firstly processed for 96 min on machine M1. It is subsequently machined for 36 min on machine M4, processed for 36 min on machine M8, and finally, it is machined for 72 min on machine M9).

**Figure 4.** Sequences of machines for each job/part.

The alternative designs that are depicted in Figure 5 can be divided into two groups: two-cell solutions (three design alternatives, i.e., Cell Designs 1–3) and three-cell solutions (two design alternatives, i.e., Cell Designs 4–5). The sequence of machines for each job is not violated in all the alternatives.

As can be seen in this figure, Cell Design 1 consists of two cells: 11 machines are located in the first cell, and 12 machines are located in the second cell. Parts P1–P4 and P7–P11 are processed in the first cell, and the rest of the parts marked as P12–P19 are processed only in the second cell. Parts P5 and P6 are partially processed in the second cell and finalized in the first cell.

Cell Design 2 contains 25 machines. There are 16 machines in the first cell and 9 machines in the second cell. Parts P1–P11 and P18 are machined in the first cell, and the rest of the parts are produced in the second cell.

Cell Design 3 is quite similar to Cell Design 2, except that machine M1 is eliminated in the second cell, and parts P15 and P16 are first produced in the first cell and subsequently finished in the second cell.

Cell Design 4 includes three cells with 8 machines in the first cell, 10 machines in the second cell, and 8 machines in the third cell. Parts marked as P1, P3, P7–P9, and P11 are machined in the first cell, while P3 is completed in the second cell. Parts P2, P4–P6, P10, P15, P16, and P18 are produced in the second cell, but P15 and P16 are finished in the third cell. Parts P12–P14, P17, and P19 are machined in the third cell.

**Figure 5.** Cell design alternatives and their material flows.

Wϭϳ

Cell Design 5 is also divided into three cells. The first cell consists of 5 machines, 15 machines are located in the second cell, and 4 machines are placed in the third cell. P1, P3, P7–P9, and P11 are machined in the first cell, while P1, P3, and P7–P9 are completed in the second cell. P2, P4–6, P10, P14–P16, and P18 are produced in the second cell, while P15 and P16 are finalized in the third cell. P12, P13, P17, and P19 are machined only in the last cell.

#### *4.2. Application of the Performance Indicators on Cell Designs*

In this sub-section, the above-mentioned indicators are applied to the five CM alternatives. Obtained operational complexity values, makespans, and production line balancing rates are summarily shown in Table 2.

The obtained values are further used as input data for the purpose of multi-criteria comparison applying the AHP method.


**Table 2.** Enumerated results of all the criteria and sub-criteria.

\* Best obtained value.

#### *4.3. Assessment of Manufacturing Cell Designs Using AHP Method*

Firstly, the hierarchical structure of the AHP method is created. The overall focus is aimed at the selection of the optimal manufacturing cell design(s) from the five cell design alternatives. For this purpose, these five alternatives are compared using the criteria shown in Figure 6.

**Figure 6.** Hierarchical structure of the AHP method for selection of an optimal manufacturing cell.

Once the hierarchy is constructed, the alternatives are pairwise compared for each of the criteria based on the preferences using the scale of relative importance (see Table 3).



Note: When comparing the best alternative with the worst alternative, the difference between them is maximally 100%. Based on this difference value, the percentages 25%, 50%, and/or 75% are proposed.

The performance of each cell design (CD) with regard to each criterion is indicated by the following pairwise comparison matrices, and at the same time, it is assumed that all criteria are equal to each other (see Figure 7).

Subsequently, the priority vectors are calculated according to the additive normalization method for all the criteria shown in Table 4.


EŽƚĞ͗ƵĞƚŽƚŚĞĨĂĐƚƚŚĂƚƉƌĞĐŝƐĞůLJĐĂůĐƵůĂƚĞĚƋƵĂŶƚŝƚĂƚŝǀĞǀĂůƵĞƐǁĞƌĞƵƐĞĚĨŽƌƚŚĞƉƵƌƉŽƐĞŽĨƚŚĞƉĂŝƌǁŝƐĞĐŽŵƉĂƌŝƐŽŶ͕ŝƚŝƐŶŽƚŵĞĂŶŝŶŐĨƵůƚŽĂƐƐĞƐƐƚŚĞŵƵƐŝŶŐƚŚĞ ĐŽŶƐŝƐƚĞŶĐLJĐŽĞĨĨŝĐŝĞŶƚ͘

**Figure 7.** Obtained values of pairwise comparison matrices for Cell Designs 1–5.

**Table 4.** Obtained priority vectors values for all the criteria.


Finally, the aggregated relative priorities from Table 4 are enumerated, which allows ranking CM alternatives, as shown in Figure 8.

**Figure 8.** Final comparison of cell designs.

As can be seen from Figure 8, Cell Design 4 is considered the optimal design.

#### **5. Conclusions**

Summarily, it can be stated that both of the three-cell solutions better satisfied the determined criteria of manufacturing cell design performance than all the three two-cell solutions. It can be empirically explained by this that cells practically represent modules, and a modular manufacturing layout design is better than an integral design.

Moreover, from the obtained results of the computational experiments, it can be noted that:

According to both complexity indicators, the three-cell solutions are less complex than the two-cell solutions. The lower complexity of the three-cell designs against the two-cell designs can be comprehended in a way that the scheduling of cell designs with a higher number of cells is less complicated than in the case with a smaller number of cells. This statement comes from the fact the probability that parts are produced on given machines is higher than in the case of the CM design with a smaller number of cells;

Based on the makespan results, the three-cell solutions better satisfied the minimization of the total time needed to finish all the jobs than the two-cell solutions;

From the viewpoint of the PLB indicator, the two-cell solutions offer better balancing of machines than the three-cell solutions.

As mentioned, the case study was taken from the work of Yan and Irani [65], who compared two-cell solutions with three-cell solutions based on selected criteria such as the number of intra-cell flows and inter-cell flows, scheduling, etc., to point out the advantages and disadvantages of two-cell solutions and three-cell solutions. Our aim was to identify the optimal cell design solution(s) from the alternatives based on a multi-criteria decision-making approach, where five selected criteria were used. The main benefit of the used approach lies in the objective approach of the pairwise comparison of cell design alternatives in the decision-making process.

Related future research could be oriented to employ other criteria to assess manufacturing cells design performance in order to bring new findings for practitioners and researchers.

**Author Contributions:** Conceptualization, V.M.; methodology, V.M. and Z.S.; validation, Z.S.; formal analysis, V.M. and Z.S.; writing—original draft preparation, V.M. and Z.S.; writing—review and editing, V.M.; visualization, Z.S. and J.N.; supervision, V.M.; project administration, V.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the KEGA Project No. 025TUKE-4/2020 granted by the Ministry of Education of the Slovak Republic.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Notations**


#### **References**


### *Article* **Meta-Heuristic Technique-Based Parametric Optimization for Electrochemical Machining of Monel 400 Alloys to Investigate the Material Removal Rate and the Sludge**

**Vengatajalapathi Nagarajan 1, Ayyappan Solaiyappan 1,\*, Siva Kumar Mahalingam 2, Lenin Nagarajan 2, Sachin Salunkhe 2, Emad Abouel Nasr 3, Ragavanantham Shanmugam <sup>4</sup> and Hussein Mohammed Abdel Moneam Hussein 5,6**


**Abstract:** Electrochemical machining (ECM) is a preferred advanced machining process for machining Monel 400 alloys. During the machining, the toxic nickel hydroxides in the sludge are formed. Therefore, it becomes necessary to determine the optimum ECM process parameters that minimize the nickel presence (NP) emission in the sludge while maximizing the material removal rate (MRR). In this investigation, the predominant ECM process parameters, such as the applied voltage, flow rate, and electrolyte concentration, were controlled to study their effect on the performance measures (i.e., MRR and NP). A meta-heuristic algorithm, the grey wolf optimizer (GWO), was used for the multi-objective optimization of the process parameters for ECM, and its results were compared with the moth-flame optimization (MFO) and particle swarm optimization (PSO) algorithms. It was observed from the surface, main, and interaction plots of this experimentation that all the process variables influenced the objectives significantly. The TOPSIS algorithm was employed to convert multiple objectives into a single objective used in meta-heuristic algorithms. In the convergence plot for the MRR model, the PSO algorithm converged very quickly in 10 iterations, while GWO and MFO took 14 and 64 iterations, respectively. In the case of the NP model, the PSO tool took only 6 iterations to converge, whereas MFO and GWO took 48 and 88 iterations, respectively. However, both MFO and GWO obtained the same solutions of EC = 132.014 g/L, V = 2406 V, and FR = 2.8455 L/min with the best conflicting performances (i.e., MRR = 0.242 g/min and NP = 57.7202 PPM). Hence, it is confirmed that these metaheuristic algorithms of MFO and GWO are more suitable for finding the optimum process parameters for machining Monel 400 alloys with ECM. This work explores a greater scope for the ECM process with better machining performance.

**Keywords:** electrochemical machining (ECM); material removal rate (MRR); nickel presence (NP); grey wolf optimizer (GWO); moth-flame optimization algorithm (MFO); Monel 400 alloys

#### **1. Introduction**

For the last few decades, electrochemical machining (ECM) has been used for machining the macro- and micro-components of tool steel, carbides, superalloys, and titanium

**Citation:** Nagarajan, V.; Solaiyappan, A.; Mahalingam, S.K.; Nagarajan, L.; Salunkhe, S.; Nasr, E.A.; Shanmugam, R.; Hussein, H.M.A.M.

Meta-Heuristic Technique-Based Parametric Optimization for Electrochemical Machining of Monel 400 Alloys to Investigate the Material Removal Rate and the Sludge. *Appl. Sci.* **2022**, *12*, 2793. https://doi.org/ 10.3390/app12062793

Academic Editors: Vladimir Modrak and Zuzana Soltysova

Received: 23 January 2022 Accepted: 1 March 2022 Published: 9 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

alloys. These materials are heavily used in the automotive and aerospace industries, as well as electronics, optics, medical devices, and communications. ECM is used as a high-priority machining tool because of its specific advantages, such as minor tool wear, mechanical forceless machining, no heat-affected zones, high surface quality, low roughness, and stressfree surface products. ECM is the reverse process of electroplating and can machine any hard materials or complicated shapes [1]. In ECM, the workpiece and tool are connected with the anodic pole and the cathodic pole of the electrochemical circuit, respectively, and the electrolyte is pumped through the inter-electrode gap (IEG) between the tool and the workpiece. A high-current DC power supply is allowed to pass through this electrochemical circuit to dissolve the metal from the workpiece in the form of metal hydroxides (i.e., sludge). ECM is an efficient and low-cost machining method for Ni-based alloys [2]. Several works were carried out to study the machinability of Ni alloys by ECM [3–5]. Ni alloys are primarily used in aircrafts, power generation turbines, rocket engines, chemical processing plants, and nuclear power plants. Poor surface traits, a short insert life, high manufacturing costs, and low productivity are associated with the machining of nickel alloys.

Monel 400 alloy is one of Ni alloys that exhibits excellent characteristics of corrosive resistance, strength, and toughness in challenging conditions [6]. When machining Ni-based alloys, ECM generates a large quantity of toxic sludge in the form of Ni and chromium hydroxides aside from gaseous byproducts, acids, sulphates, nitrates, oils, and metal ions. Environmentally sustainable machining of Ni-based alloys in ECM is needed while maintaining productivity and quality [7]. However, the discharge of nickel ions into the electrolyte slurry is significant. Therefore, it is necessary to use ECM with proper cutting conditions so that it produces sludge with fewer toxic byproducts. This hazardous emission may be retained in the electrolyte by repeated use of the electrolyte. However, such a process might decrease the machining performance. Thus, it is proposed to have optimized ECM process parameters to produce a better MRR and fewer hazardous material emissions. To the best of the authors' knowledge, investigating the amount of nickel present in the sludge during electrochemical machining of Monel 400 alloys has not been attempted before. Optimizing the process parameters is an essential task in the ECM process, as the optimum process parameters improve the performance and economics of machining. The selection of ECM process parameters is carried out based on the experience and expertise of the machinist or machining handbooks. The process parameters selected based on the operator's experience rarely assure high efficiency and quality machining. The machining handbooks can be a handy choice for a few applications only. In most cases, the selected optimum process parameters are far from the best, and this hampers the best utilization of ECM. Selecting the optimum values for the process parameters without proper optimization requires elaborate, costly, time-consuming, and tedious experimentation. Therefore, researchers have used different soft computing techniques to optimize the ECM process parameters in various operations [8]. Mukherjee and Chakraborty implemented a biogeography-based optimization (BBO) algorithm in ECM to determine the optimum machining conditions [9]. The BBO algorithm was inspired by the migration behavior of species among the habitats in nature. A genetic algorithm (GA) was integrated with a desirability function (DF) for simultaneous optimization of multiple objectives of the ECM process [10]. The particle swarm optimization (PSO) algorithm was adopted as an optimization strategy for determining the optimum parameters [11]. The PSO algorithm mimics the concept of the social interaction of birds in a flock. To perform the multi-objective optimization, the cuckoo optimization algorithm (COA) was used with the objective function, which was derived from the adaptive neuro-fuzzy inference system (ANFIS) model [12]. The COA was developed based on the obligate brood parasitism of some cuckoo species. A multi-objective Jaya (MO-Jaya) algorithm was implemented in the ECM process to obtain multiple optimal solutions [13]. The differential evolution (DE) algorithm has been used for both the single- and multi-objective optimization of process parameters [14]. Singh and Shukla (2017) considered the ECM process to evaluate the performance of the black hole algorithm (BHA). The performance measures of the MRR and overcut (OC) were taken

to obtain the optimum machining parameters using the BHA [15]. The harmony search (HS) algorithm was combined with the DF so that an HS-DF optimizer was proposed to optimize the ECM process parameters [16]. Diyaley and Chakraborty conducted a comparative analysis with meta-heuristics such as the firefly algorithm (FA), DE algorithm, ant colony optimization (ACO) algorithm, and teaching-learning-based optimization (TLBO) algorithm for the optimization of various control parameters of an ECM process [17]. Apart from these meta-heuristics, artificial neural networks (ANNs) and fuzzy logic (FL) were also employed to model the ECM experiments and subsequent optimization process [18,19]. Several optimization algorithms have been reported in the literature, and their efficiency was tested with many benchmark functions and real-time applications. The grey wolf optimizer (GWO), a recently developed meta-heuristic, was proposed by Mirjalili et.al [20]. The GWO is inspired by the hunting behavior of the grey wolves in nature, and it ensures a better trade-off between the exploration and exploitation abilities of the algorithm. It was applied in many engineering applications, emphasizing its efficiency [21]. The significant characteristic of this algorithm is that it does not converge toward some local optima and helps store the best possible solutions obtained so far by its social hierarchy nature [22]. Kharwar and Verma implemented GWO to minimize milling performances (i.e., MRR, cutting force, and surface roughness) and proved the application potential in a manufacturing environment [23]. Omkar and Shalaka optimized the wire electrical discharge machining (WEDM) parameters with the GWO algorithm [24]. Shankar and Ankan experimented with the efficiency of the GWO algorithm for parametric optimization of the abrasive water jet machining (AWJM) process and found it to be successful [25].

A performance comparison of the response surface methodology (RSM), genetic algorithm (GA), and GWO algorithm was conducted for prediction of the surface roughness in ball-end milling of hardened steel, and GWO was found to be superior [26]. The grey relational analysis (GRA)-based GWO was implemented for milling experiments to determine the machining conditions of the spindle speed, feed rate, and depth of cut for the optimum surface roughness, cutting force, and MRR [27]. Mirjalili proposed another meta-heuristic, the moth-flame optimization (MFO) algorithm, which was inspired by the navigation method of moths in nature [28]. It gained attention for solving engineering optimization problems immediately. The improved MFO algorithm was proposed by Li et al. to maintain a balance between global and local searching [29]. A milling optimization problem was solved to emphasize its effectiveness in solving manufacturing problems [30]. The MFO algorithm was implemented to identify the optimal set of turning parameters to minimize the machinability indices. Its effectiveness was compared against the genetic algorithm (GA), grasshopper algorithm, GWO, and PSO algorithms and to be found superior to the other algorithms [31]. The optimization of the surface roughness in deep hole drilling was performed with the MFO algorithm [32]. The MFO algorithm was employed to optimize the process parameters of plasma arc cutting (PAC) of Inconel 718 superalloy [33]. Based on the literature on GWO, this algorithm is not utilized much in manufacturing optimization problems. These reviews conclude that the recently developed intelligent techniques are more advanced, and their scope of implementation in broader manufacturing applications can be explored. It is also noted that the efficiency of these new evolutionary algorithms (i.e., GWO and MFO algorithms) were not explored in the ECM process. The famous PSO algorithm is still being used to optimize machining parameters [34,35]. Hence, PSO is considered a benchmark algorithm to compare GWO and MFO algorithms in terms of convergence, computational time, and the actual Pareto optimal front.

#### **2. Objectives**

To determine the optimum ECM parameters for machining Monel 400 alloys for the maximum MRR and minimum nickel toxic emission in the electrolyte, meta-heuristics such as GWO, MFO, and PSO algorithms were employed.

#### **3. Materials and Experimentation**

The commercially available nickel-based alloy Monel 400 was used as a test specimen in this investigation. The chemical composition of the Monel 400 alloys is shown in Table 1 [36]. The machinability of Monel 400 alloy is very difficult, as it work hardens during machining. Therefore, it is a tough task to machine these alloys using conventional machining techniques.



The experiments were carried out on the eco-friendly electrochemical machining (EECM) tool. The schematic diagram of the EECM tool is shown in Figure 1. EECM comprises a power supply system, electrolyte supply system, filtering system, tool feed mechanism, work holding and position system, control panel, frame, and housing. The electrolyte (e.g., NaCl aqueous solution) is allowed to flow at the rate between 1 L/min and 10 L/min through the inter-electrode gap (IEG) of 0.1–0.6 mm. The direct current (DC) potential (11–15 V) with a current density of 20–110 A/cm2 is applied across IEG between a cathodic copper tool and an anodic workpiece. At low current densities, the MRR is low [37]. The higher current density of 110 A/cm2 was used to attain high current efficiency. According to Faraday's laws of electrolysis, the anodic dissolution rate of workpiece depends on the electrochemical properties of the workpiece metal, the electrolyte properties, and the power supply conditions. The tool is fed against the anodic workpiece, which is firmly fixed on the vice. The machining conditions of IEG, voltage (V), and current (C) are set in the control panel. The filtration system consists of layers of bio-absorbents compartments and membrane filter. A photograph of the indigenously developed EECM is shown in Figure 2. The tool's diameter and electrolyte exit hole were 8 mm and 2 mm, respectively. The operating conditions of EECM are summarized in Table 2.

**Figure 1.** Schematic diagram of EECM.

**Figure 2.** EECM experimental set-up.



#### *Experimental Design and Measurements*

The MRR is directly proportionate to the feed rate of the tool [36]. A high feed rate leads to electrolyte boiling or choking in the tool and workpiece gap [37]. Therefore, the tool feed rate was set to the possible minimum of 0.1 mm/min to have stable movement of the tool through the workpiece. When the applied voltage is high, the current machining increases to a high MRR. A low voltage results in poor machining performance [10,11]. The voltage range of 11–16 V was considered in many such investigations in ECM to have better performance in terms of the MRR and surface finish [38,39]. Therefore, the voltage range (10–15 V) available in the set-up was used for this investigation. When the voltage is raised above a particular level, it increases the hydrogen gas bubble generation at the tool electrode, resulting in high resistivity of the electrolyte and a decrease in the current density at the workpiece [40–43]. A high concentration decreases the mobility of ions, which results in poor anodic dissolution. The concentration levels in the range of around 100–200 g/L influence the MRR significantly [10]. According to the literature, the main process parameters governing the ECM are the voltage (V), flow rate, and EC. Therefore, these parameters were considered in this investigation, and their levels are shown in Table 3.


**Table 3.** Process parameters and levels.

The IEG was set to 0.1 mm, and a current of 50 A was set in the DC rectifier. The machining was performed for 5 min. The central composite design (CCD) of the response surface methodology (RSM) was adopted with the help of Minitab-17 software for the process parameters, and its levels are shown in Table 3. In this work, 20 experimental tests were carried out with different parameter combinations as shown in Figure 3. According to the concept of CCD, the center point combination was repeated 5 more times during experimentation to reduce the pure error.

**Figure 3.** Central composite design for the experimentation.

A HCl concentration of 1% was added to the NaCl electrolyte to minimize the sludge formation in the IEG. The electrolyte was pneumatically pumped from a stainless steel reservoir. A 191E water analysis kit (Environmental & Scientific Instruments Co, Haryana, India.) was used to observe the electrolyte's pH, conductance, and temperature.

The sludge discharged from ECM was tested with atomic absorption spectroscopy (Agilent Technologies, Bangalore, India), which is shown in Figure 4. The results show the nickel content of 250 mg/L in the sludge. The precipitate of the electrolyte was tested via energy dispersive X-ray spectroscopy (EDAX, JEOL-JSM-6390, Chennai, India) and confirmed the presence of nickel. Figure 5 represents the EDAX results of the electrolyte sludge. The peaks in the figure confirm the presence of more nickel, chloride, and sodium particles. Therefore, in this invetigation, the experimental results for MRR and Nickel Presence (NP) were observed and tabulated in Table 4.

**Figure 4.** Atomic absorption spectroscopy.

**Figure 5.** EDAX spectrum for nickel in precipitate.


**Table 4.** Experimental design and performance results.

#### **4. Optimization Techniques and Procedures**

The following details describe the development of mathematical models representing EECM characteristics for the MRR and the amount of nickel present, as well as the GWO and MFO techniques for optimizing the process parameters.

#### *4.1. Mathematical Modeling of Experimentation*

To study the effect of the EECM process parameters on the MRR and the discharge of nickel elements into the electrolyte, a regression model was developed to calculate the response value (as the output) in terms of different parameters (as the input). In this work, the quadratic form of the regression equation, shown in Equation (1), was established using the experimental data. The regression coefficients obtained for both the MRR and NP are shown in Table 5:

$$\begin{array}{c} \text{CV}\_{ij} = a\_j + b\_j \ \* \mathcal{X}\_{1i} + c\_j \ \* \mathcal{X}\_{2i} + d\_j \ \* \mathcal{X}\_{3i} + e\_j \ \* \mathcal{X}\_{1i} \* \mathcal{X}\_{2i} + f\_j \ \* \mathcal{X}\_{1i} \* \mathcal{X}\_{3i} + g\_j\\ \ast \mathcal{X}\_{1i} \* \mathcal{X}\_{3i} + h\_j \ \* \mathcal{X}\_{1i} ^2 + i\_j \ \* \mathcal{X}\_{2i} ^2 + j\_j \ \* \mathcal{X}\_{3i} ^2 \end{array} \tag{1}$$

where *CVij* is the predicted value of the *i*th trial for the *j*th response, *a, b,* ... *, j* are the coefficients of the variables, and *X1, X2*, and *X3* are independent variables.



The regression models for the given objectives can be formulated as in Equations (2) and (3):

$$\begin{array}{l} MRR = -2.381 & + 0.0123 \ast EC + 0.112 \ast V + 0.621 \ast FR + 0.0003 \ast EC \ast V\\ & - 0.003 \ast EC \ast FR + 0.0143 \ast V \ast FR - 0.00003 \ast EC^2\\ & - 0.007 \ast V^2 - 0.065 \ast FR^2 \end{array} \tag{2}$$

$$\begin{array}{c} NP = & 389 - 0.551 \ast EC - 41.8 \ast V - 49.5 \ast FR + 0.0011 \ast EC \ast V + 0.131 \\ \ast EC \ast FR + 1.23 \ast V \ast FR + 0.00154 \ast EC^2 + 1.561 \ast V^2 \\ \ast 4.64 \ast FR^2 \end{array} \tag{3}$$

#### *4.2. TOPSIS Method for Multi-Objective Optimization*

The Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS) is a method used for converting multiple objectives into a single objective, which is consequently applied in optimization techniques for decision making [44,45]. TOPSIS is based on choosing the alternative with the shortest geometric distance from the positive ideal solution and the longest geometric distance from the negative ideal solution. By implementing the TOPSIS method, multiple objectives are converted into a single objective value. Algorithm 1 shows the the pseudo-code for TOPSIS method.

**Algorithm 1** Pseudo-code for TOPSIS.

1: Read alternate and objectives matrix—*Oij* with weights (*Wj*) and types of objectives (*otj*)

2: For each alternative (i = 1, 2, 3, . . . , m) and objective (j = 1, 2, 3, . . . , n)

3: Compute Normalized value of *Oij* using *Nij* <sup>=</sup> *Oij* ∑*<sup>m</sup> <sup>i</sup>*=<sup>1</sup> *O*<sup>2</sup> *ij* 4: Calculate Performance matrix (*Aij*) using *Aij* = *Nij* ∗ *Wj* 5: End 6: For each objective (j = 1, 2, 3, . . . , n)

7: Determine positive ideal (*Pj*) and negative ideal solution (*Mj*) 

8: For minimization objective—*Pj* = min 1<*i*<*m Aij* and *Mj* = max 1<*i*<*m Aij* 9: For maximization objective—*Pj* = max *Aij* and *Mj* = min *Aij*

1<*i*<*m*

10: End

11: For each alternative (i = 1, 2, 3, . . . , m)

12: Calculated Ideal (*SPi*) and negative ideal separation (*SMi*)

$$SP\_{l} = \sqrt{\Sigma\_{j-1}^{n} \left(A\_{l\bar{j}} - P\_{\bar{l}}\right)^{2}} \\ SM\_{l} = \sqrt{\Sigma\_{j-1}^{n} \left(A\_{l\bar{j}} - M\_{\bar{l}}\right)^{2}}$$

13: Determine relative closeness (*Ri* ) 14: End

15: Rank alternatives w.r.t. *Ri* in descending order

#### *4.3. Grey Wolf Optimizer*

The grey wolf optimizer (GWO) was inspired by the hunting behavior of grey wolves in nature. Grey wolves usually live in a pack of 5–12 members with their dominant social hierarchy [20]. The alpha (α) wolf, the pack leader, makes decisions about hunting. The beta (β) wolf helps the alpha wolf in making decisions. The omega (ω) wolves, the lowest in ranking, are responsible for submitting information to all the wolves. Others, called delta (δ)

1<*i*<*m*

wolves, should respect alpha and beta wolves and dominate the omega wolves. The three steps in the hunting process of grey wolves are (1) tracking the prey, (2) encircling the prey, and (3) attacking the prey. GWO is simulated with the mathematical representation of the hunting process of the grey wolves to solve complex engineering problems. The optimal solution to the problem is considered the prey. To mathematically model the hunting process of these grey wolves when designing GWO, the alpha (α) wolf is considered the fittest solution. Consequently, the second- and third-best solutions are represented by beta (β) and delta (δ) wolves, respectively. The remaining candidate solutions are considered to be omega (ω) wolves.

The first step in the hunting process is encircling the prey. The mathematical formulation to mimic the encircling process for the prey and wolf are represented in the Equations (4) and (5):

$$
\stackrel{\rightarrow}{D} = \left| \stackrel{\rightarrow}{\mathcal{C}} \* \stackrel{\rightarrow}{X}\_p(t) - \stackrel{\rightarrow}{X}(t) \right| \tag{4}
$$

$$
\overrightarrow{X}\ (t+1) = \overrightarrow{X}\_p\ (t) - \overrightarrow{A} \* \overrightarrow{D} \tag{5}
$$

where <sup>→</sup> *<sup>D</sup>* is a distance vector between the position of the prey and grey wolf. <sup>→</sup> *<sup>C</sup>* and <sup>→</sup> *A* are coefficient vectors, *<sup>t</sup>* is the current iteration, and <sup>→</sup> *Xp* and <sup>→</sup> *X* are the positions of prey and a randomly chosen grey wolf, respectively [20]. The coefficients <sup>→</sup> *<sup>C</sup>* and <sup>→</sup> *A* are determined using the following formulas:

$$
\stackrel{\rightarrow}{\vec{A}} = 2\stackrel{\rightarrow}{\vec{a}} \* \stackrel{\rightarrow}{r\_1} - \stackrel{\rightarrow}{\vec{a}} \tag{6}
$$

$$
\overrightarrow{\dot{C}} = 2 \ast \overrightarrow{r\_2} \tag{7}
$$

where <sup>→</sup> *a* is a decrease from 2 to 0 linearly over the course of iterations and <sup>→</sup> *<sup>r</sup>*<sup>1</sup> and <sup>→</sup> *r*<sup>2</sup> are random vectors (0, 1).

A grey wolf is in the position of <sup>→</sup> *<sup>X</sup>*, updating its position <sup>→</sup> *X* (*t* + 1) according to the position of the prey <sup>→</sup> *Xp*. The position can be updated with respect to the current position by adjusting the values of the <sup>→</sup> *<sup>A</sup>* and <sup>→</sup> *C* vectors. This process enables the GWO to search the n-dimensional solution space of the given problem efficiently. The alpha (α), beta (β), and delta (δ) wolves guide the hunting process, and their positions are represented mathematically as in Equations (8)–(14), while the other wolves (ω) update their positions randomly:

$$
\stackrel{\rightarrow}{D}\_{\infty}^{\rightarrow} = \left| \stackrel{\rightarrow}{\mathcal{C}}\_{1} \* \stackrel{\rightarrow}{X}\_{\infty} - \stackrel{\rightarrow}{X} \right| \tag{8}
$$

$$
\stackrel{\rightarrow}{D}\_{\emptyset}^{\rightarrow} = \left| \stackrel{\rightarrow}{\mathcal{C}}\_{1} \* \stackrel{\rightarrow}{X}\_{\emptyset} - \stackrel{\rightarrow}{X} \right| \tag{9}
$$

$$
\stackrel{\rightarrow}{D}\_{\mathcal{S}} = \left| \stackrel{\rightarrow}{\mathcal{C}}\_{1} \* \stackrel{\rightarrow}{X}\_{\mathcal{S}} - \stackrel{\rightarrow}{X} \right| \tag{10}
$$

$$
\stackrel{\rightarrow}{X}\_1 = \stackrel{\rightarrow}{X}\_{\infty} - \stackrel{\rightarrow}{A}\_1 \* \stackrel{\rightarrow}{D}\_{\infty} \tag{11}
$$

$$
\stackrel{\rightarrow}{X}\_2 = \stackrel{\rightarrow}{X\_{\beta}} - \stackrel{\rightarrow}{A\_2} \* \stackrel{\rightarrow}{D\_{\beta}} \tag{12}
$$

$$
\stackrel{\rightarrow}{X}\_{3} = \stackrel{\rightarrow}{X}\_{\delta} - \stackrel{\rightarrow}{A}\_{3} \* \stackrel{\rightarrow}{D}\_{\delta} \tag{13}
$$

$$
\stackrel{\rightarrow}{X}(t+1) = \frac{\stackrel{\rightarrow}{X}\_1 + \stackrel{\rightarrow}{X}\_2 + \stackrel{\rightarrow}{X}\_3}{3} \tag{14}
$$

As the value of decreases, the search radius of the grey wolves reduces, and they get closer to the prey and attack over the last iterations of the GWO. This ensures the proper trade-off between exploration and exploitation by focusing on exploration initially and exploitation in the last iterations.

#### *4.4. Moth-Flame Optimization Algorithm (MFO)*

The moth-flame optimization (MFO) algorithm takes inspiration from the navigational mechanism (i.e., transverse orientation) of moths during night flights. The moth flies by maintaining a fixed angle relative to the moon. Since the moon is so far away from the moth, this mechanism assures the linear flight of moths. Although a lateral orientation is practical, moths are often trapped to repeatedly circle many close artificial or natural point light sources until they are exhausted. This behavior of moths is shown in Figure 6.

**Figure 6.** Spiral flight path of moths near artificial or natural point light sources [20].

Mathematical Formulation of the MFO Algorithm

According to the MFO algorithm, the moth is considered the candidate solution to the problem. The variables are represented by the positions of the moths in the solution space. Moths can fly in the solution space by changing the position vectors. Since this algorithm is a swarm-based intelligence optimization algorithm, the position vectors of the moth population can be represented in the following matrix (Equation (15)):

$$\mathbf{M} = \begin{bmatrix} m\_{1,1} & \cdots & m\_{1,d} \\ \vdots & \ddots & \vdots \\ m\_{n,1} & \cdots & m\_{n,d} \end{bmatrix} \tag{15}$$

where *n* is the number of moths and *d* represents the number of dimension variables to be solved. The following array represents the list of fitness value vectors corresponding to all moths:

$$\text{OM} = \begin{array}{c} \text{OM}\_1 \\ \vdots \\ \text{OM}\_n \end{array} \tag{16}$$

Each moth updates its position with the corresponding unique flame to avoid the algorithm being trapped in the optimal local value, which significantly supports the algorithm's exploration ability. Hence, the flame positions in the solution space corresponding to its moths are represented in the following array (Equation (6)):

$$\mathbf{F} = \begin{bmatrix} f\_{1,1} & \cdots & f\_{1,d} \\ \vdots & \ddots & \vdots \\ f\_{n,1} & \cdots & f\_{n,d} \end{bmatrix} \tag{17}$$

For these flames, the corresponding column of the fitness value vectors is represented as follows:

$$\text{OF} = \begin{vmatrix} OF\_1 \\ \vdots \\ OF\_n \end{vmatrix} \tag{18}$$

The position update strategies for both the moth and flame matrices are different during the iteration process. The moths are the solution points that move within the solution space, and the flames are the best solution points obtained by moths iteratively so far. The best solutions of the moths are updated to the positions of the flames in the next iteration. With this approach, the MFO algorithm can find the optimal global solution. The position update behavior of each moth relative to a flame is expressed mathematically by the following equation:

$$M\_{\hat{\imath}} = S\left(M\_{\hat{\imath}\prime}F\_{\hat{\jmath}}\right) \tag{19}$$

where *Mi* is the *i*th moth, *Fj* is the *j*th flame, and S is the spiral function.

The spiral function's initial point should start from the position of the moth, and the endpoint should be at the position of the flame. The range fluctuation of the spiral function should not exceed the boundaries of the search space. According to the above conditions, the logarithmic spiral function of the moth's flight path is defined as in Equation (20):

$$S(M\_i, F\_{\bar{j}}) = D\_i \ast e^{bt} \ast \cos(2\pi \text{t}) + F\_{\bar{j}} \tag{20}$$

where *Di* is the linear distance of the *i*th moth for the *j*th flame and *b* is an index of the shape of the logarithmic spiral. The magnitude of *t* is represented by Equation (21), where *a* is represented by Equation (22):

$$t = (a - 1) \* rand + 1\tag{21}$$

The path coefficient *t* ∈ [*r*, 1] represents the distance between the moth and the flame's position in the next optimization iteration. The variable *r* decreases linearly from −1 to −2 as the number of iterations in the optimization iteration increases. The coefficient from *t* = −1 to −2 represents that the moths' position is close to the flame, and *t* = 1 shows that the moths are farther from the flame:

$$a = -1 + I \\ iteration\* \left(-\frac{1}{T\_{\max}}\right) \tag{22}$$

*Di* is expressed as follows:

$$D\_i = \left| F\_{\bar{j}} - \mathcal{M}\_{\bar{i}} \right| \tag{23}$$

The above equations ensure the algorithm's balanced global and local search capabilities. When the value of *t* is smaller, the moth converges to the flame. As the moth gets closer to the flame, its position around the flame is updated more quickly. After each iteration, the flames are updated based on fitness values. In the next iteration, the moth updates its position according to the updated sequence of the flame. The first moth updates its position with the flame's best fitness value, and the last moth updates its position with the worst fitness value in the list. Hence, the flame position matrix *F* contains the optimal solutions of the current iteration. The number of flames can be reduced adaptively in the iterative process to balance the algorithm's global search and local search capability in the solution space as follows:

$$frame\_{no} = round\left(N - l\*\frac{N-1}{T}\right) \tag{24}$$

where *l* is the current iteration number, *N* is the initial maximum number of flames, and *T* represents the maximum number of iterations.

#### **5. ANOVA and Parametric Influence on Performances**

Analysis of variance (ANOVA) on the process parameters and the effects of the main process parameters on the performance measures was performed. This and Fisher's test (F-ratio) were performed to justify the adequacy of the developed regression models as in Equations (2) and (3). Values of "Prob > F" less than 0.05 indicate the model's significance and terms. Values greater than 0.1 indicate that the model terms are insignificant.

The model F-values in Table 6 indicate that the developed models were very significant. The two-way interaction terms such as EC\*V, EC\*EC, and V\*FR for the MRR and NP models were insignificant when the values of "Prob > F" were greater than 0.05. The terms V\*V and EC\*FR for the MRR and NP models, respectively, were not significant, but the coefficient of determination values for both models were good, as shown in Table 7. The normal probability curves in Figure 7 confirm the adequacy of the model equations. The high coefficient of determination values (R2 = 90.34% for the MRR and R2 = 92.29% for the NP) for both responses confirm its adequacy for optimization study.

**Table 6.** Analysis of variance for quadratic models of MRR and NP.


**Figure 7.** Normal probability plots for MRR and NP.


**Table 7.** Model summary.

#### *Parametric Influence on the Performance Measures*

Figure 8 depicts the surface plots of the MRR and NP responses for electrochemical machining of Monel 400 alloys. In each surface plot, one parameter was kept constant, and other parameters were varied with defined levels. The surface plots show that the machining parameters V, FR, and EC influenced the material removal and the discharge of nickel elements in the sludge significantly. It can be observed that the change in electrolyte concentration from 150 g/L to 190 g/L gradually increased the sludge removal. A higher concentration helped increase the speed of the electrochemical reactions and resulted in more removal (Zhang, 2010; Ayyappan and Sivakumar, 2015). A lower concentration led to a lower ion content, which reduced the rate of the electrochemical reaction (Trimmer, Hudson, Kock, and Schuster, 2003). A flow rate above 1.5 L/min removed the sludge from the gap between the tool and the workpiece better and exposed the specimen's new surface for different electrochemical reactions (Ayyappan and Sivakumar, 2014). High flow rates increased the nickel ions in the electrolyte discharge and the MRR (Bhattacharyya and Sorkhel, 1999).

Figure 9 represents the main effect plots of the MRR and NP responses using MiniTab17 software. It was noticed that the electrolyte flow rate and electrolyte concentration were the influencing factors for the MRR and NP, respectively. Variations in the voltage also produced significant results in the responses.

Figure 10 shows the interaction plots of the ECM parameters for different process attributes. As shown, the interaction terms such as EC\*V and V\*FR had a lean influence on the responses. The MRR and NP models were insignificant when the "Prob > F" values were more significant than 0.05. The term EC\*FR influenced the NP at a high flow rate of 3 L/min.

**Figure 8.** *Cont*.

**Figure 8.** Surface plot of responses for various combinations of parameters.

**Figure 9.** Main effects plots of MRR and NP.

**Figure 10.** Interaction plots for MRR and NP.

#### **6. Multi-Objective Optimization of ECM Process Parameters**

*6.1. RSM Optimization Tool*

Since the RSM as a DOE strategy was employed in this work, the optimum process parameters were obtained with its built-in optimization tool. The RSM optimization tool used the steepest ascent (SA) method combined with the desirability function (DF) for finding the optimum parameters. The variables and response constraints were set as shown in Table 8.

**Table 8.** Limits of process and performance parameters.

**Table 9.** Multiple response prediction results using RSM.


The optimization results are shown in Table 9, and the optimization plot is shown in Figure 11.


**Figure 11.** Response optimization plot.

#### *6.2. Applications of Meta-Heuristics for Optimization of the Process Parameters*

To apply the various meta-heuristic algorithms such as MFO, GWO, and PSO in this work, the TOPSIS method was adopted to convert the multiple objectives into a single objective. A performance comparison of these algorithms was carried out. The Algorithms 2–4 show the pseudo-codes of the MFO, GWO, and PSO respectively.

**Algorithm 2** Grey Wolf Optimization Algorithm


**Algorithm 3** Moth-Flame Optimization Algorithm

1: Initialize the parameters for Moth-flame 2: Initialize Moth position *Mi* randomly 3: For each *i = 1:n* Calculate the fitness valute *fi* 4: End 5: While *(i* ≤ *imax)* 6: Update the position of *Mi* 7: Calculate the no. of flames 8: Compute the fitness value *fi* 9: If *(i = 1)* then *F = sort (M) OF = sort (OM)* 10: Else *F = sort (Mt-1, Mt) OF = sort (Mt-1, Mt)* 11: End 12: For each *i = 1:n* 13: For each *j = 1:d* Update the values of *r* and *t* 14: Calculate the value of *D* w.r.t. corresponding Moth 15: Update *M(i,j)* w.r.t. corresponding Moth 16: End 17: End 18: End 19: Print the best solution


The terms in GWO and their equivalent meanings in the optimization problem are listed in Table 10.


**Table 10.** Terms in GWO and optimization problems.

Similarly, the equivalent terms for the MFO and PSO algorithms for an optimization problem are shown in Table 11. The algorithm parameters are shown in Table 12.




**Table 12.** Algorithm parameters for GWO, MFO, and PSO algorithms.

The steps in implementing the GWO algorithm in this experimentation are shown in Figure 12. The flowchart for implementing the MFO algorithm in the present work is shown in Figure 13.

**Figure 12.** Implementation of GWO algorithm.

**Figure 13.** Implementation of MFO algorithm.

The PSO implementation is shown in Figure 14. The efficiency of these algorithms is compared in the convergence plot of Figure 15. In maximization of the MRR, the MFO convergence plot was attained after the 60th iteration, while GWO and PSO performed better. In minimization of the NP, PSO converged well before MFO and GWO.

The optimum process parameters obtained by different meta-heuristics and the RSM are tabulated in Table 13. The convergence plot of the Pareto optimal solutions is shown in Figure 16, and the Pareto optimal set of process parameters out of 12 runs is presented in Table 14.

**Table 13.** Optimum process parameters and the response values.


**Figure 14.** Implementation of PSO algorithm.

**Figure 15.** Comparison of convergence plots for various algorithms for MRR and NP.

**Figure 16.** Pareto optimal solutions for various algorithms.


**Table 14.** The Pareto optimal solutions and their responses for various algorithms.

Both the MFO and GWO algorithms produced the better trade-off performances with the same solutions of EC = 132.014 g/L, V = 13. 2406 V, and FR = 2.8455 L/min for MRR = 0.242 g/min and NP = 57.7202 PPM.

#### **7. Conclusions**

In the present study, ECM operations were carried out on Monel 400 alloys by varying the applied voltage (11–15 V), flow rate (1–3 L/min), and electrolyte concentration (130–190 g/L). The meta-heuristic algorithms (i.e., moth-flame optimization (MFO) algorithm, grey wolf optimizer (GWO), and particle swarm optimization (PSO)) were implemented to find out the optimum process parameters for producing the best performance measures for the material removal rate (MRR) and nickel presence (NP) in the sludge. The regression model equations were developed to determine the optimum parametric combination. The high coefficient of determination values (R2 = 90.34% for MRR and R2 = 92.29% for NP) for both responses confirmed their adequacy. The algorithms were tuned to obtain the best possible feasible solution. This study proved that the MFO and GWO algorithms could be successfully utilized to find the best process parameter setting for the ECM process for Monel 400 alloys. The TOPSIS algorithm was implemented to convert multiple objectives into a single objective used in meta-heuristic algorithms. It was observed from the surface, main, and interaction plots of this ECM experimentation that all the process variables significantly affected the performances. The PSO algorithm converged very quickly in 10 iterations, while GWO and MFO took 14 and 64 iterations, respectively, for convergence of the MRR. In the case of the NP model, the PSO tool took only six iterations to converge, whereas MFO and GWO took 48 and 88 iterations, respectively. However, both MFO and GWO obtained the better trade-off performances with the same solutions of EC = 132.014 g/L, V = 13. 2406 V, and FR = 2.8455 L/min for MRR = 0.242 g/min and NP = 57.7202 PPM. The Pareto optimal set of solutions provided by the GWO, MFO, and PSO algorithms for the ECM process was beneficial to the industries, as it contained a wide range of optimal values. A particular solution from the Pareto optimal set can be chosen according to the specific needs of the objectives. Hence, it was confirmed that the metaheuristic algorithms of MFO and GWO were suitable for finding the optimum process parameters for machining Monel 400 alloys with ECM. Future studies may be explored by combining the artificial neural network (ANN) and fuzzy logic (FL) concepts with MFO and GWO to determine the optimum ECM process parameters.

**Author Contributions:** Conceptualization and methodology, V.N., A.S., L.N. and S.K.M.; writing original draft preparation, V.N., A.S., L.N. and S.K.M.; writing—review and editing, S.S., E.A.N., R.S. and H.M.A.M.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors extend their appreciation to King Saud University for funding this work through Researchers Supporting Project number (RSP-2021/164), King Saud University, Riyadh, Saudi Arabia.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Optimization of Process Parameters for Turning Hastelloy X under Different Machining Environments Using Evolutionary Algorithms: A Comparative Study**

**Vinothkumar Sivalingam 1, Jie Sun 1, Siva Kumar Mahalingam 2, Lenin Nagarajan 2,\*, Yuvaraj Natarajan 2, Sachin Salunkhe 2, Emad Abouel Nasr 3, J. Paulo Davim <sup>4</sup> and Hussein Mohammed Abdel Moneam Hussein 5,6**


**Abstract:** In this research work, the machinability of turning Hastelloy X with a PVD Ti-Al-N coated insert tool in dry, wet, and cryogenic machining environments is investigated. The machinability indices namely cutting force (CF), surface roughness (SR), and cutting temperature (CT) are studied for the different set of input process parameters such as cutting speed, feed rate, and machining environment, through the experiments conducted as per L27 orthogonal array. Minitab 17 is used to create quadratic Multiple Linear Regression Models (MLRM) based on the association between turning parameters and machineability indices. The Moth-Flame Optimization (MFO) algorithm is proposed in this work to identify the optimal set of turning parameters through the MLRM models, in view of minimizing the machinability indices. Three case studies by considering individual machinability indices, a combination of dual indices, and a combination of all three indices, are performed. The suggested MFO algorithm's effectiveness is evaluated in comparison to the findings of Genetic, Grass-Hooper, Grey-Wolf, and Particle Swarm Optimization algorithms. From the results, it is identified that the MFO algorithm outperformed the others. In addition, a confirmation experiment is conducted to verify the results of the MFO algorithm's optimal combination of turning parameters.

**Keywords:** Hastelloy X; turning; cutting force; surface roughness; liquid nitrogen; grass-hooper optimization algorithm; moth-flame optimization algorithm

#### **1. Introduction**

Nickel-based (Ni) alloys attract more researchers nowadays for their broader applications in the fields like aerospace, automobile, biomedical, and allied industries. Hastelloy is one of the Ni-based alloys, and it holds few unique characteristics like good strength-toweight ratio, resistance to corrosion, higher melting temperature, good toughness, etc. [1]. Mainly, Hastelloy X is used to fabricate the combustion chamber of an aircraft engine because of its high heat-resisting property. However, the holding of all the above-said

**Citation:** Sivalingam, V.; Sun, J.; Mahalingam, S.K.; Nagarajan, L.; Natarajan, Y.; Salunkhe, S.; Nasr, E.A.; Davim, J.P.; Hussein, H.M.A.M. Optimization of Process Parameters for Turning Hastelloy X under Different Machining Environments Using Evolutionary Algorithms: A Comparative Study. *Appl. Sci.* **2021**, *11*, 9725. https://doi.org/10.3390/ app11209725

Academic Editor: Vladimir Modrak

Received: 6 September 2021 Accepted: 11 October 2021 Published: 18 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

properties by Hastelloy X, resulting in very poor machinability. In this sense, the manufacturing industries face a difficult task in improving Hastelloy X machinability using traditional machining methods. [2]. Furthermore, the reduction of cutting forces (CF), surface roughness (SR), and cutting temperature (CT) during Hastelloy X machining adds to the difficulty of achieving good machinability. As a result, several researchers have worked on various research projects over time to increase the machinability of Hastelloy. Furthermore, they performed these tests under dry, wet, and cryogenic cooling conditions in order to demonstrate an increase in machinability. Therefore, these literatures are critically reviewed, and the extracted information is given here for ready reference to the readers.

Kadirgama et al. [3] studied the impact on cutting force by the parameters, namely axial depth, cutting speed, and feed rate while milling Hastelloy C-22HS. The models using Response Surface Methodology were developed using experimentation and Finite Element Analysis to predict the optimized cutting force. Kadirgama et al. [4] investigated the tool behavior such as tool wear and tool life during machining of Hastelloy C-22HS under wet conditions. PVD and CVD multilayer coated carbide tools were used for machining. The tool life was decreased in all the cases while increasing the cutting parameters, namely cutting speed (*v*c), feed rate (*f*), and axial depth (ap). Altin [2] studied the machinability of Ni-based (Hastelloy X) alloy under dry cutting conditions. The CF and SR were analyzed against the multilayer coated insert and various *v*c. The experimentation results showed that the abrasiveness of the carbide particles on the tool and the mechanical loading had a growing influence on the CF. Sofuo ˘glu et al. [5] studied the impact of the *v*c, tool extended length, and novel methods, namely Conventional Turning (CT), Ultrasonic Assisted Turning (UAT) and Hot-Ultrasonic Assisted Turning (HUAT) on the SR, ap, and CT while machining Hastelloy X. The reduction in SR and increment in regular ap and CT were attained in UAT and HUAT compared to CT. Dhananchezian [6] conducted the machinability study on Hastelloy C-276 under dry and cryogenic liquid nitrogen (LN2) cooling conditions using turning operation. The output responses such as CT, CF, SR, chip morphology, and tool wear under dry turning were compared with LN2 cooling-based turning. A considerable reduction in all the output responses was noted under liquid nitrogen cooling-based turning.

Kesavan et al. [7] conducted the CNC turning of Hastelloy C276 by varying *v*<sup>c</sup> and the fixed values of *f*, ap. The experimentation was executed under dry and LN2 conditions. Further, Deform 3D analytical tool was used to create the simulation model based on the experimental design to identify the optimal cutting conditions. From the experimentation and simulated model results, it was evident that the cutting temperate and machining forces have been significantly reduced while machining under cryogenic cooling conditions rather than dry conditions. Dhananchezian and Rajkumar [8] examined the SR and Tool Wear characteristics of Nimonic 90 alloy and Hastelloy C-276 dry turning. During the turning process, various cutting inserts were used. In both cases, the roughness and tool wear metrics were observed to be larger as the turning length was increased. Dhananchezian and Rajkumar [8] made a comparative analysis on the tool wear rate and SR during the turning of Hastelloy C-22 underneath dry and LN2 cooling conditions. A substantial drop in the SR was found in the turning of Hastelloy under LN2 cooling rather than dry turning.

Oschelski et al. [9] used the Box-Behnken method to design the experiments by considering the ranges of parameters, namely *v*c, ap, lubricating conditions, constant *f*, and (wet, dry, and reduced quantity lubrication) for finish turning the Hastelloy X. The experimental results showed that the *v*c, ap, and interactions were the most significant factors affecting the SR. Next, Venkatesan et al. [10] reported the machinability study on Hastelloy X with PVD and CVD coated tools in comparison with dry and Minimum Quantity Lubrication (MQL) conditions. A mixture of coconut oil with Hexagonal Boron Nitride (HBN) nanoparticles was used as nanofluid for lubrication. Significant reductions in CF, SR, and tool wear were observed in MQL-PVD combination than MQL-CVD and dry-PVD. Finally, Sivalingam et al. [11] investigated the influence of whisker-reinforced

ceramic tools on tool wear, SR, and tool chattering under dry and Atomization-based Cutting Fluid (ACF) cooling conditions when turning Inconel 718 material. Investigation results stated that the flank wear and SR of the tool were significantly reduced under ACF cooling conditions due to limited notching and fracture of the tool edge at the tool-chip interface.

Zhao et al. [12] investigated the characteristics of chip formation when machining NiTi shape memory alloys under different *v*c with constant *f*, ap. The shape of the chip and microstructure were examined to expose the chip flow behavior. The martensitic phase transformation seemed to have a noticeable effect on the material flow behavior and indeed on the chip formation. [11] investigated the possibility of improving the machinability of Inconel 718 alloy under a dry and atomized spray cutting fluid system. The turning of Inconel 718 alloy with ceramic inserts was carried out by varying the cutting parameters. The output responses such as tool wear, power consumption, surface topography, machine vibrations, chip morphology, and machining cost were analyzed against the experimental design of input parameters. It was observed that the atomized spray cutting fluid technique yielded better results than dry machining.

The effect of LN2 cooling in improving the machinability of Hastelloy X is discussed in the following literature. Chetan et al. [13] investigated the turning of Nimonic 90 alloy using uncoated tungsten carbide inserts under the modes like dry, MQL, and cryogenic cutting. At lower *v*c, the cutting performance of the cryogenically treated tool was good than the untreated tool. But, the performance of the tool under MQL and LN2 was good in terms of minimum tool wear at a higher cutting speed. Further, a good SR was obtained under dry and MQL modes than LN2 cooling mode at all levels of cutting speed. Iturbe et al. [14] compared the effects of liquid nitrogen and MQL based cryogenic cooling with conventional cooling. For short machining times, the cryogenic cum MQL cooling outperformed conventional cooling.

Sivaiah and Chakradhar [15] compared the results of LN2 machining like tool wear, feed force, CF and CT, chip characteristics, and SR with the wet condition during machining of heat-treated 17-4 Precipitation Hardenable Stainless Steel. The LN2 machining outperformed even at high *f* to reduce all the above-said parameters compared with wet machining. Tebaldo et al. [16] studied the machinability of Inconel 718 under different machining conditions and lubricating systems. The highest wear resistance was obtained while using the CVD-coated tools under conventional lubricated conditions. But, the MQL system provided good lubrication than cooling with lesser cost and low environmental impact. Shokrani et al. [17] investigated the impact of using different cooling systems, namely MQL, cryogenic and hybrid of cryogenic and MQL, during the CNC milling of Inconel 718 alloy material. Comparatively, the hybrid cooling system yielded better results in terms of good machinability, less SR, and greater tool life. Mehta et al. [18] studied the parameters such as SR, CF, and tool wear during machining of Inconel 718 material. During machining, various sustainable environments, namely dry state, MQL, LN2 cooling, hybridization of cold air and MQL, and hybridization of MQL and LN2, were used. The input parameters such as ap, *f*, and *v*c were kept constant during machining under all the above-said environments. Better surface finish and minimum cutting force were observed during the cold air and MQL environment. Alternatively, the very least tool wear was observed under MQL and LN2 hybrid cutting environment than the dry environment.

Further, the researchers had used different optimization tools to identify the suitable process parameter values for minimizing the manufacturer's objectives. A few of them are discussed here. Khalilpourazari and Khalilpourazary [19] proposed an algorithm, namely Robust Grey Wolf Optimizer (RGWO), to minimize total production time by identifying the optimal input parameters multi-pass milling process. The parameter tuning during optimization was carried out using the Taguchi method. Further, an efficient constraint handling approach was implemented to handle the complex constraints of the problem. The results concluded that the RGWO outperformed the meta-heuristic algorithms such as the multi-verse optimizer and dragonfly algorithm and the other solution methods

in the literature. Khalilpourazari and Khalilpourazary [20] developed the lexicographic weighted Tchebycheff method to obtain the optimal decision parameters of the grinding process for maximizing the quality of the surface and production rate and minimizing the machining time and cost. GAMS software was used for this purpose. Khalilpourazari and Khalilpourazary [21] used a novel strategy, namely Robust Stochastic Novel Search, to identify the optimal values of the grinding process parameters to minimize process cost and to maximize the rate of production and surface quality. The proposed method outperformed previously proposed methodologies and novel algorithms, including Multi-Population Ensemble Differential Evolution and Heterogeneous Comprehensive Learning Particle Swarm Optimization.

Rao et al. [22] optimized the abrasive water-jet machining parameters to minimize the kerf and surface roughness using the Jaya algorithm and the multi-objective Jaya algorithm. Better results were obtained through the used algorithms than the simulated annealing, particle swarm optimization, firefly algorithm, cuckoo search algorithm, blackhole algorithm, and biogeography-based optimization algorithms. Further, the PROMETHEE method was used to handpick a specific solution among the possible Pareto-optimal solutions obtained through the proposed algorithms based on the given requirements. Rao et al. [23] obtained the optimal set of process parameters of focused ion beam micro-milling, laser cutting, wire-electric discharge machining, and electrochemical machining processes. The maximization of material removal rate and minimization of SR were considered as the objectives in all the processes. The multi-objective Jaya algorithm was implemented to find the optimal solutions in all the cases. The test results showed that the implemented algorithms produced good results compared to other algorithms such as Genetic Algorithm, Non-dominated Sorting Genetic Algorithm, iterative search, and biogeography-based optimization algorithm, Khalilpourazari and Khalilpourazary [24] carried out the optimization of grinding process parameters to improve the SR and reduce production cost and time. A multi-objective dragonfly algorithm was employed for optimizing the process parameters. Results revealed that the proposed algorithm outperformed the Non-dominated Sorting Genetic Algorithm-II. Khalilpourazari and Khalilpourazary [25] proposed a novel hybrid algorithm, Sine–Cosine Whale Optimization Algorithm, to optimize the process parameters of the multi-pass milling process by minimizing the total production time. Almeida et al. [26] optimized the variable-angle composite cylinders via filament winding manufacturing process using GA. Similarly, Wang et al. [27] proposed a reliability-based design optimization technique to improve the buckling load of winding cylinders subjected to radial compression. The moving search windows in the Kriging metamodel are used to accelerate its convergence and reduce the number of training iterations. The results of this study demonstrated the advantages of adopting a variable stiffness design for achieving a maximum load capacity. Almeida et al. [28] proposed a genetic algorithm (GA) to enhance the strength of a cylindrical shell under internal pressure by optimizing the stacking sequence. The results offered asymmetric and non-conventional angles for internally pressured composite tubes, as opposed to the well-known ± 55◦ winding angle advice (for first ply failure approach).

In this research work, turning experiments are conducted on the Hastelloy X material using the PVD TiAlN carbide insert tool under dry, wet, and LN2 environments. The *v*c, *f*, ap, and machining environment are considered input turning process parameters, and CF, SR, and CT are considered machinability indices. The evolutionary algorithms namely grasshopper optimization (GHO) [29–31], genetic algorithm (GA) [32–34], particle swarm optimization (PSO) [35,36], moth flame (MFO) [37,38], grey wolf optimization (GWO) [39–41] algorithms are used to identify the optimal set of turning process parameters (MATLAB R2020b version). A clear picture of the experimentation and subsequent processes are detailed in Figure 1.

**Figure 1.** Flow chart of experimental part and metaheuristic algorithm.

#### **2. Experimentation Details**

The Hastelloy X bar having 20 mm diameter and 300 mm length is used for conducting the turning experiments on a C6140H turning machine. VNMG160408-SM1105 PVD TiAlN cutting tool inserts are used to do the turning operation. The impact of tool wear on the machinability indices is completely eliminated by using new inserts every time. The CF (Tangential Force, Fz) is calculated with a 9257B Kistler dynamometer, and the value is manipulated with dynoware software. The workpiece surface roughness (Ra) is measured using a contact-type surface roughness tester (TR200), a cutoff length of 0.8 mm, and a traverse length of 4 mm. For measuring the CT, a FORTIC 226 infrared imaging sensor is used. Cryogenic equipment consists of a self-pressurized pump, cryogenic dewar tank capacity of 50 L. At a pressure of 0.3 bar, LN2 was sprayed onto the work-tool interface using a copper nozzle diameter of 3 mm. Figure 2 depicts a schematic representation of the experimental setup. The detailed experimental conditions are shown in Table 1.

**Figure 2.** (**a**) Schematic View, Experimental Setup (**b**) Dry (**c**) Wet (**d**) Cryogenic Machining.


**Table 1.** Experimental Conditions.



In this work, turning experiments were performed using L<sup>3</sup> <sup>27</sup> full factorial experimental design using Minitab 17. three factors are considered for this experiment: *v*c, *f*, and environment; each factor has three different levels, as shown in Table 2.



The experimental design and the corresponding measurement of machinability indices are presented in Table 3.


**Table 3.** Experimental design values.


**Table 3.** *Cont.*

#### **3. Results and Discussion**

This section is divided into three case studies. Case study 1: Minimization of machinability indices individually; Case study 2: Simultaneous minimization of dual machinability indices by considering three combinations; Case study 3: Simultaneous minimization of all three indices. The quadratic Multiple Linear Regression Models (MLRM) are formulated for evaluating the minimum values of machinability indices in all the cases. The Moth-Flame Optimization (MFO) algorithm is proposed to identify the optimal set of turning process parameters in view of minimizing the objectives. The effectiveness of the proposed algorithm is tested against the results of other optimization algorithms such as Genetic Algorithm, Grass-Hooper Optimization (GHO), Grey-Wolf Optimization (GWO), and Particle Swarm Optimization (PSO). Pseudocode for optimization algorithms is shown in the Figure 3. The general parameters used in algorithms are the maximum population size: 50, and the maximum no. of iterations: 100 (MATLAB R2020b). Twenty-seven runs are executed for each algorithm in all the cases. The evaluated results from the case studies are discussed below.

#### *3.1. Case Study 1*

Cutting force analysis and its minimization plays a crucial role in machining operation and understanding the cutting phenomena of the work material in different environments (dry, wet, and LN2 machining). Moreover, after machining, the work material must be superior in surface quality [33]. On the other hand, minimizing the machining zone's temperature is necessary to retain the cutting temperature as low as possible. In this case study, the MLRM for minimizing all the machinability indices is developed using Minitab 17. The developed MLRM are given in Equations (1)–(3).

**Objective functions 1.**

$$\begin{array}{l}\text{Minimize Cutting} = 212.7 - 0.244A + 703B + 20\text{C} \\ \text{Force} & -0.00571A^2 + 6289B^2 - 8.11\text{C}^2 \\ & -0.60AB + 0.053AC - 143BC \\ R^2 = 0.98, \text{ }Adj \, R^2 = 0.98 \end{array} \tag{1}$$

**Objective functions 2.**

$$\begin{array}{l}\text{Minimize } Surface = 3.56 - 0.008A - 9.65B + 0.73C\\\text{Roughness} & + 0.000015A^2 + 19.3B^2 - 0.265C^2\\ & + 0.023AB + 0.00028AC - 0.65BC\\\ R^2 = 0.98, \ Add R^2 = 0.97\end{array} \tag{2}$$

#### **Objective functions 3.**

$$\begin{array}{l}\text{Minimize } \text{Cutting} = -0.083 + 0.076A + 2521B + 341.21 \text{C} \\ \text{Temperature} \qquad + 0.0033A^2 - 1400B^2 - 118 \text{C}^2 \\ \qquad - 1.42AB - 0.16AC - 396.67BC \\ \text{ $R^2 = 0.97, \text{ $ Adj $ R^2 = 0.96}$ } \end{array} \tag{3}$$

Three responses were considered in this section R1, R2 and R3, CF, SR and CT, respectively, as shown in Table 4.

**Table 4.** Minimization of machinability indices using evolutionary algorithms for Case study 1.


Figure 4a present the convergence plot for cutting force using the evolutionary algorithms. It is observed that the response is converged in iteration no. 2 using the MFO algorithm. Provided the same response value is converged in iteration no. 3, 10, 61, and 78 based on GA, PSO, GHO, and GWO algorithms, respectively. Further, the response value is very minimum in the MFO algorithm and maximum in GA. The simplicity of the MFO algorithm along with the speed in searching is the prime reason for obtaining the best results in the present work. It is additionally inferred that the CF is very minimum in the LN2 environment compared to the dry and wet machining environments. The LN2 nozzle spray droplet acts as cushioning effect of the machining zone and, thus, minimizes the *v*c [42].


**Figure 3.** Pseudocode for optimization algorithms (**a**) MFO (**b**) GHO (**c**) GA (**d**) GWO (**e**) PSO.

**Figure 4.** Convergence plot for case study 1 (**a**) Cutting Force (**b**) Surface Roughness (**c**) Cutting Temperature.

Similarly, the convergence plots for the other two machinability indices are given in Figure 4b,c, respectively. The order of preferences of algorithms based on their performance in obtaining the minimum response values is MFO-GWO-GHO-GA-PSO for SR and MFO-GHO-PSO-GWO-GA for CT. Chattered vibrations are generally existing in the machining process, and this is significantly affecting the SR. The SR values are lower in the cryogenic machining than the dry and wet machining from the experimental value and predicted data.

Further, the cutting temperature is directly related to CF and SR. CT increases with an increase in *v*c. The flow of cryogenic LN2 (−196 ◦C) between tool and workpiece interface greatly reduces the CT in the machining zone compared with dry and wet machining [42].

#### *3.2. Case Study 2*

In this study, the simultaneous minimization of dual machinability indices with three different combinations using evolutionary algorithms is considered. The Pareto front analyses of all the algorithms with respect to CF vs. SR, CT vs. CF, and CT vs. SR are given in Figure 5a–c, respectively. Further, the TOPSIS method is used to convert the dual machinability indices into a single objective. Hence, the global minimum value of the machinability indices is obtained using TOPSIS results (Figure 5d–f) for all the evolutionary algorithms. In addition to that, the performance of evolutionary algorithms is validated using the hypervolume indicator. From these Figure 5, it is inferred that the MFO algorithm outperformed others in all three cases. The results are presented in Table 5.

**Table 5.** Minimization of machinability indices based on pareto front analysis and TOPSIS method for Case study 2.


The general observations are when CF increases, and SR worsens after machining; When CF increases, TE increases significantly and affects the tool life; and when TE increases, then SR decreases. Further, the LN2 environment minimizes the CF, CT, and SR due to a fine droplet of LN2 acting as a film barrier in the tool and workpiece interface, thus reducing chatter vibration and built-up edge formation. This is used to improve the tool life too.

**Figure 5.** Pareto optimality (**a**–**c**) Multi-objective for single run (**d**–**f**) Multi-Objective for 27 runs.

#### *3.3. Case Study 3*

In this case study, the simultaneous minimization of all three machinability indices is carried out using evolutionary algorithms. As in case study 2, the TOPSIS method is used to convert all three machinability indices to a single objective. Further, the quality indicators, namely Diversity (DIV), Inverted Generational Distance (IGD), and Hyper Volume (HV), are used to opt out the best evolutionary algorithm which provides minimum values of machinability indices along with the corresponding set of turning process parameters. The algorithm, which has a higher DIV and HV value and a lower IGD value, is to be considered as the best algorithm among others. The DIV and IGD values are calculated using Equations (4) and (5), respectively. The hypervolume is calculated based on the Pareto analysis. The calculated values of quality indicators for all the algorithms are presented in Table 6.

*DIV* = ∑*k <sup>j</sup>*=<sup>1</sup> (*f* max *<sup>j</sup>* −*<sup>f</sup>* min *<sup>j</sup>* )<sup>2</sup> (4)

$$IGD = \frac{\sqrt{\sum\_{i=1}^{n} d\_i^2}}{n} \tag{5}$$

$$d\_{\vec{i}} = \sqrt{\sum\_{j=1}^{no} \left(o\_{\vec{i}\vec{j}} - o\_{b\vec{j}}\right)^2} \tag{6}$$

where,

*Oij*—*i*th run *j*th objective value; *Obj*—Best *j*th objective value; *di*—Euclidean distance.

**Table 6.** Minimization of machinability indices based on IGD, DIV and HV for Case study 3.


Further the statistical analyses of the quality indicators are carried out using Minitab software to study the performance and consistency of the algorithms. The normal probability plots for the quality indicators IGD, DIV, and HV are given in Figure 6a–c, respectively, and the summary reports of the same are given in Figure 6d–f, respectively. From Figure 6d–f and as well as the findings (higher values of DIV and HV and Lower value of IGD) from Table 5, it is concluded that the MFO algorithm outperformed others.

The reasons for the better performance of MFO compared with other algorithms are explained here. The GA became gradually a dominant optimization technique compared to deterministic approaches mainly due to the higher probability of local solutions avoidance. However, the main drawback of GA was the stochastic nature of this algorithm which resulted in finding different solutions in every run. Despite the relatively high convergence rate of PSO, it has the drawback of premature convergence to local optimal and ineffectiveness in exploring the whole search space. Similarly, the original version of the GWO algorithm has the drawbacks of low solving accuracy, bad local searching ability, and slow convergence rate. Similarly, the disadvantage of GHO was being easy to fall into the local optimum which prevented the search process from finding a better solution. On the other hand, MFO is able to locate the local and global optimal solutions accurately with less computational time [43].

**Figure 6.** (**a**–**c**) Normal Probability plot (**d**–**f**) summary report.

#### **4. Conclusions**

The purpose of this study was to minimize machinability indices CF, SR, and CT while performing the turning of Hastelloy X. Three levels of turning process parameters namely cutting speed (*v*c), feed rate *f*, and machining environment were considered for performing the experiments under L27 orthogonal array basis. Further, the MFO algorithm was used to identify the optimal set of turning process parameters to minimize the machinability indices individually and simultaneously. Three case studies were carried out for this purpose. The conclusion drawn from these case studies is given below.


Based on the results of all three case studies, the MFO algorithm effectively predicted the optimal set of turning process parameters in view of minimizing the machinability indices individually and simultaneously when compared with other algorithms. Further, the other machinability indices such as tool life and machining cost will also be considered in addition to the existing indices as the future work.

**Author Contributions:** Conceptualization, V.S.; Methodology, V.S., J.S. and Y.N.; Experimental design, V.S. and J.S.; Experimental setup, V.S.; Measurements, V.S., S.K.M. and Y.N.; Investigation, L.N., S.K.M. and Y.N.; Resources, L.N., S.K.M. and Y.N.; Visualization, V.S., L.N., S.K.M. and S.S.; Writing—Original Draft Preparation, V.S., J.S., L.N., S.K.M., Y.N., S.S., E.A.N., J.P.D. and H.M.A.M.H.; Writing—Review & Editing, V.S., J.S., L.N., S.K.M., Y.N. and S.S.; Supervision, J.S. and V.S.; Project administration, V.S. and J.S.; Funding acquisition, J.S. and S.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has received funding from King Saud University through Researchers Supporting Project number (RSP-2021/164), King Saud University, Riyadh, Saudi Arabia. Moreover, this is part of the project was financially supported by Fundamental Research Fund of Shandong University under the funding code No. 2019HW040. the Future for Young Scholars of Shandong University, China under the funding code No. 31360082064026.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors extend their appreciation to King Saud University for funding this work through Researchers Supporting Project number (RSP-2021/164), King Saud University, Riyadh, Saudi Arabia.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Reflections on the Customer Decision-Making Process in the Digital Insurance Platforms: An Empirical Study of the Baltic Market**

**Gedas Baranauskas \* and Agota Giedre Raišien ˙ e \* ˙**

Institute of Leadership and Strategic Management, Mykolas Romeris University, Ateities Str. 20, LT-08303 Vilnius, Lithuania

**\*** Correspondence: gedasbaranauskas@mruni.eu (G.B.); agotar@mruni.eu (A.G.R.); Tel.: +370-62-151-887 (G.B.)

**Abstract:** Multifold effects of the COVID-19 global health crisis and economic lockdowns are reflected in the insurance industry, and are predicted to expand to the post-COVID-19 era. It is expected that, within a short period of time, the current worldwide situation, in regards to the coronavirus pandemic, will be reflected in new trends, regarding customer behavior, organizational management, and culture, as well as reveal improved business management models, legacy infrastructure, and service systems in insurance organizations. Here, a focus on end-user preferences, data, and their behavior modeling in digital platforms are major practical drivers within the modern insurance concept, but there is a paucity of researches within the theoretical synthesis of consumer decisionmaking (CDM) models, information system theories, and behavioral economics concerning modern insurance-specific value chains and digitalized decision-making processes. Therefore, the following research aims to expand upon the existing scientific knowledge of end-user behavioral patterns and process frameworks in the Baltic insurance market, by including and examining a factor group of technological enablers and a digital environment. Research results in digitalization, personalization, and customization levels within the Baltic non-life insurance market are homogenous with a leading position of Estonia and overall evaluations ranging between Satisfied and Rather Good. There are also three major factor groups and process stages identified, which influence insurance purchase decision-making in digital insurance platforms in the Baltic market.

**Keywords:** digital platforms; decision-making; insurance; Baltic; customization; personalization

#### **1. Introduction**

Countries across the world underwent an unprecedented global health crisis and subsequent lockdowns in 2020–2021 due to the coronavirus disease (COVID-19) pandemic, moreover, social and economic outcomes are still hard to predict at a full scope. It is argued that insurance organizations should focus on the development of hybrid servicebased and customer-driven business models, sustainable and innovative digital products, and consolidations, with InsurTech and technology companies, for better digitization of distribution channels [1,2]. The ongoing global pandemic has revealed a lack of flexibility, emotional connection, and data harmonization with the day-to-day needs of end-users, as well as personalized, situation-based, and easily customizable insurance services, as a part of practical improvements toward the modern insurance [3,4]. It is important to note intensive and dynamic changes in the market structure, which are mostly determined by a high competition among traditional peers and new virtual incumbents, non-traditional insurance service providers, such as "Big Tech", and product manufacturers [5,6].

From a theoretical point of view, the research refers to predecessor studies in the field of non-life insurance consumer decision-making, such as Hsee and Kunreuther (2000) [7], Kunreuther and Pauly (2006; 2015) [8,9], and a cross-border study report by the European Commission (2017) [10]. The present research also follows key findings and limitations

**Citation:** Baranauskas, G.; Raišiene,˙ A.G. Reflections on the Customer Decision-Making Process in the Digital Insurance Platforms: An Empirical Study of the Baltic Market. *Appl. Sci.* **2021**, *11*, 8524. https:// doi.org/10.3390/app11188524

Academic Editors: Vladimir Modrak and Zuzana Soltysova

Received: 4 August 2021 Accepted: 6 September 2021 Published: 14 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of scientific studies on behavioral patterns in the Baltic insurance market, such as research by Kiyak and Pranckeviˇciut¯ e (2014) [ ˙ 11], and longitudinal studies conducted by Ulbinaite et al. ˙ (2010, 2011, 2013) [12–15]. Additionally, novel studies of comparison websites [2], multi-sided platforms [1] in the insurance industry and the Service-Dominant Architecture (SDA) [16], as well as insurance literacy [16,17], were considered. However, a recent increase in practical demands for hyper-personalized support services and digitally customized insurance outcomes of on-demand and usage-based insurance reveal a research gap. There is a need for a multi-tier factor evaluation of the digital insurance environment and multi-agent-based model simulations oriented toward decision making in digital insurance platforms. Therefore, this research aims to extend the existing scientific knowledge on insurance end-user behavioral patterns and frameworks by including and examining factors of customization, personalization, and technological enablers. To investigate this situation, the following research questions have been raised:


From a methodological point of view, the research follows a triangulation logic, combining an online survey with Fuzzy and Likert logic, simplified art-based method outcomes for information presentation, and an embedded, descriptive case study with statistical data analysis methods, which summarize key findings of the online survey. The structure of this research paper is composed of five main sections. The first section introduces the theoretical background of the topic. The second section provides a detailed description of the research and survey methodology and applied methods. Results of the online survey are presented and analyzed in the third section. The fourth section presents a discussion on research limitations and possible future research directions. The fifth section concludes the main findings from theoretical and practical analyses.

#### **2. The Theoretical Background of Digital Non-Life Insurance Concepts and Insurance Customer Decision-Making Processes**

Although the importance of non-life insurance products for the financial wellbeing of individuals and society has been recognized for a long time, a practical spread is still considered vague, i.e., the lack of insurance literacy, financial decision-making skills, and a number of underinsured persons are observed [17]. On a theoretical level, it is argued that modern insurance-related decisions and customer behavior cannot be sufficiently explained by applying traditional neoclassical economic and financial theories and ignoring cognitive and emotional factors.

In general, fundamental studies on consumer decision-making in insurance services, insurers, and market behavior, were carried out by Hsee and Kunreuther in 2000 [7] and Kunreuther and Pauly in 2006 [8]. However, the phenomenon of insurance purchases under these studies has been analyzed from a holistic standpoint, resulting in a focus on high-level benchmark models of insurance demand and supply, legal regulation, and government insurance market principles. The interdisciplinary nature of the behavior of consumers in the insurance market was analyzed, mostly regarding monetary financial risks and positive theory perspectives, with limited attention to the influence of behavioral finance and economics theories, as well as heuristic and bias phenomena. This research gap was only partially covered in the last decade via a comprehensive analysis from the European Commission (2014) [10] study, and Kunreuther and Pauly in 2015 [9]. Here, assumptions on behavioral finance and economy and cognitive psychology domains have contributed to the foundation of the modern insurance concept by introducing the influence of individual risk perceptions, heuristic reflections in framing effects, and intuitive and deliberative thinking along general concepts of probability and risk, distribution channels, and socialdemographic differences in Europe [9,10,17]. Nevertheless, afterwards, an increase in scientific attention to the phenomenon of insurance digitization and digitalization has

been noticed, as the majority of research was mostly oriented to countries from Western and Central European regions, while analyses of the Baltic region (Lithuania, Latvia, and Estonia) resulted in low volumes, and were mostly limited to a one country-level (Lithuania or Latvia). A significant part of the existing scientific studies on digital insurance have been completed at a high-generality level, scattered to different parts of the insurance-specific value chain or technological solutions, without sufficiently considering the influence of modern mass customization and personalization concepts and marketing theories. A similar situation can be identified when analyzing consumer decision-making in insurance services, which, in the case of the Baltic region, are geographically limited, practically outdated, and, as in the case of Lithuania, the majority of studies were completed from 2008 to 2014, without considering customization, personalization, or digitalization factors. Therefore, the multidimensional and cross-border studies of digital insurance decisionmaking, insurance literacy, and platform business models in the Baltic region are required both theoretically and practically.

To conclude, the modern insurance combines derivatives from several concepts and theories:


Modern insurance specific value chain and co-creation activities depend on features of the core offering-purchase stage, but are also shaped by numerous internal and external factors: customer sourcing processes and an organization role within customer operations, high-level factors of sociodemographic and socioeconomic characteristics, insurance literacy, platform design and information layout, and social media and reality [12,13,17,19,24,25]. Figure 1 summarizes key elements and processes of insurance decision-making.

**Figure 1.** Features and background of the insurance decision-making process. Source: composed by the authors, by following [5,12,13,17].

Figure 1 illustrates the main process stages and influential driving factors, which affect the insurance decision-making and customer transition through the stages. In general, all influential driving factors can be divided into two groups. The first group is formed by internal factors, which are related to personal evaluation of an insurance need, financial affordability, and past and ongoing experience with insurance service providers. The second group compounds external factors, which define an objective and holistic evaluation of insurance decision-making from legal, marketing, and technical points of view. The financial and marketing-oriented factors have key role in the pre-purchase stage, where the final decision to purchase is made. It should also be outlined that technological factors and the digital environment have a significant influence on the modern insurance value chain and customer experience management. Contrary to the legal, informational, relational, or social-demographic driver factors, they are recognized equally important to the entire value co-creation process in all three stages of the insurance decision-making process. An intensive development of the internet technology, mobile devices, and digital platform business model have become critical enablers to balance the quality of service delivery, insurance personalization, and customization, with different customer expectations and experience levels:


This theoretical setup reflects the modern non-life insurance meaning, which stands for a continuous (but not a simultaneous) sequence of the insurance-related processes and multiple internal and external influencing factors within the decision-making. The complex and non-linear non-life insurance market nature requires a decentralized, digitalized, and individual-centric approach to organizational management, product configuration, and customer experience management [12,13,17]. Otherwise, digital technologies, technological advancements, and global social networks increased the recognition of non-life insurance products and services, and have been rapidly embedded into daily operations and a value co-creation chain of insurance. These trends indicate the shift of traditional insurance concept boundaries and the non-life insurance market from a traditional, provider-centric management approach and service blueprint model, where a homogenous focus on the operational efficiency and cost-service level dichotomy is dominant [26,27]. The emerging prominence of the customer-centric philosophy and holistic design framework and availability of individual-level real-time data factors reflect in newly developed, innovative, and personalized touchpoints, on customer values and experience management [26,27]. Practically, it reflects on advanced online self-service platforms, which work in combination with techniques of personalization and customization and create a new cognitive framework. It embraces an intrinsic security need for customers and leads to a simpler creation, usage, and exchange of insurance knowledge and information [5,26,28].

It is also agreed that the last two decades of the globalization process, intensive digitization, and digitalization, and rapid socioeconomic changes, have brought about a remarkable, two-fold impact on the financial service market, organizations, and customers, including the insurance industry [29]. Therefore, it is recognized that existing integrated and sophisticated quantitative analysis techniques, financial and economic theories, and consumer decision-making (CDM) models are not fully sufficient to be applied in the modern, digitalized financial service, oriented to the platform business model [30]. The analysis and modeling of the optimal financial decision-making process should be supported by the application of behavioral reasoning theory (BRT), forecasting, and multicriteria decision approach and methods with qualitative data, behavioral, and cognitive factor evaluations [29,31]. New theoretic models and conceptual frameworks should reflect both technical specifications and capabilities of digital platforms as well as recent behavioral patterns of fully digital end-users, which is no longer a linear progression through process stages, but more an iterative decision-making process [30]. The traditional CDM and hybrid decision-making models (HDCM) have been widely applied in practice since 1960s, but there is a paucity of research within the theoretical synthesis of CDM and HDCM models, concerning a modern insurance-specific value chain and digitalized decision-making process. A brief overview of the main CDM models and their influence points on digitalized insurance decision-making processes is presented in Table 1.


**Table 1.** The influence of the main CDM models on digitalized insurance decision-making processes. Source: created by the authors, following [30,32–36].


**Table 1.** *Cont.*

For a long period of time, main CDM models defined in Table 1, together with the expected utility theory have been accepted as a dominant paradigm by economic and marketing researchers, and used practically by organizations to determine and model customer responses to products or services, purchase offerings, and motivational appeals [33,36]. On the other hand, modeling of the financial attitude and decision-making process in digital platforms has specific circumstances to consider. The stage of the initial problem, problem recognition, or need arousal, together with a non-standard combination of high psychological and monetary risks, internal and external biases associated with certain consequences of financial decision-making, have a strong influence in this field [30,34]. These particular characteristics of financial decision-making require reframing processes and variables of existing CDM and HCDM models. It is also essential to identify and apply the fundamentals of customer behavior models to the comprehensive examination of the insurance service customer's engagement and characteristics. The following essential points of traditional consumer decision-making models can remarkably contribute to a conceptual framework of the digital insurance decision-making process. Fundamental outcomes of Nicosia's (1966) [37] model, as a separation of the buying process to the multiple stages, iterative and constant connection between the organization and the customer (in the form of a feedback area and repurchase cycle) can be identified in the insurance. Relevant factors, considering various endogenous and exogenous variables, marketing stimuli components, process options of partiality, and postponement of decision-making, were presented in the EKB (1968) [38] and J. A. Howard and J. N. Sheth (1969) [39] models. These elements support the content and the context of the modern insurance and decision-making process, but are also limited to apply fully, due to the logic of the linear process and the unmeasured relationships between variables [36,46]. Later revisions and elaborations of these limitations in the choice set model of S. Um and J. L. Crompton (1990) [40] and the hybrid choice model of J. Walker and M. Ben-Akiva (2002) [43] introduce relevant factors of possible disturbances, irrational behavior, and an unreliable memory of customers, as well as the initial stage of the awareness set [32,33,36]. Furthermore, simplified models by P. Kotler (1997) [43], P. Kotler and K. L. Keller (2006, 2012) [44,45], and J. E. McCarthy, W. E. Perreault, and P. G. Quester (1997) [42] reflect early complex models of buyer behavior, and in this way, support the holistic approach to the decision-making process. Overall, these models outline the pre-purchase and post-purchase stages and customer-orientated activities of identification differentiation and transformation of target groups, which are also important parts of the modern insurance concept.

#### **3. Materials and Methods**

The present research is a continuation and a supplementary portion of the authors' multi-step research (2020, 2021) within a practical status, regarding insurance digitalization, customization and personalization level, and digital insurance platforms in the Baltic non-life insurance market. The research also contributes to the research field, providing empirical validation of a conceptual framework of the digital insurance customer decisionmaking process. Therefore, data collection and analyses were carried out by using a convergent parallel research design and the following triangulation of scientific methods:

1. A descriptive thematic analysis and synthesis of scientific findings within digital non-life insurance and decision-making models.


It is argued that a convergent parallel research design and triangulation of scientific methods allow simultaneous collection, analysis, and interpretation of quantitative data and qualitative evaluations in a single research study. Moreover, complementary data can be used and transformed to a more holistic understanding of the research phenomenon [50,51]. The details of the questionnaire are presented in Table 2.

**Table 2.** Presentation of the questionnaire structure, content, and methodological foundation. Source: composed by the authors.


The questionnaire was formulated by using a combined full-blown Likert scale and a Fuzzy set of 9 points with numerical and linguistic values due to several methodological reasons:


• Provision of responses grounded with a well-balanced and gradual assessment for further interpretation within methods of descriptive statistics, and comparative and correlation analysis [54].

The selection of the fundamental scale of nine points (from 1 to 9) follows the multicriteria decision making (MCDM) approach and the logic of the analytic hierarchy process (AHP). It determines a comparison of criteria and their weight importance within the suggested conceptual framework [55]. It is also argued that a combined full-blown Likert and fundamental AHP 1 to 9 scale is relevant under the cross of the threshold at n = 10 criteria [56].

The online survey was conducted through multiple forms, distribution channels, and periods. The detailed summary of the surveying process is as follows:


The research sample was formulated by following a logic of probability sampling and a simple sample type. The probability sampling selection was related to the specific research and target population, a methodological set of selection criteria to participate in the survey, and the overall possibility of reducing the sample bias and gathering higher, unbiased data quality. Authors accept key limitations and risks of using simple random sampling—a high level of the standard error of estimate [57]. The target population was defined by using two critical selection criteria: (a) working positions and experiences in the non-life insurance field, and (b) the workplace as an insurance service provider, physically located and operating in the Baltic region (Lithuania, Latvia, or Estonia).

#### **4. Results**

A combination of statistical techniques, involving descriptive statistics and factor analyses, are recommended within the applied research to simplify an interpretation of quantitative variables and, overall, to ensure a comprehensive examination of the empirical dataset [58]. This type of multi-level qualitative and quantitative analysis allows researchers to examine the validity of theoretical constructs and the consistency of research instruments, sampling adequacy, and underlining structures and relationships among latent variables (factors) [58,59]. The selected forms of the descriptive statistics involved a manual calculation of a simple and grouped frequency distribution and measures of central tendency. The factor analysis was conducted with the SPPS statistical software (version 26) by completing exploratory factor analysis (EFA), confirmatory factor analysis (CFA), and Pearson's correlation analysis. The logical sequence of statistical analysis procedures is presented in Figure 2.

**Figure 2.** Logical sequence and content of statistical analysis. Source: created by the authors.

#### *4.1. Descriptive Statistics*

Application of the descriptive statistics provides a valuable summary of the dataset characteristics and research sample in a manageable and organized format [60]. It also contributes to the research conclusions and validation of the conceptual framework by revealing the influence of sociodemographic factors. The main characteristics of the research sample are presented in Table 3.


**Table 3.** Characteristics of the research sample. Source: composed by the authors.

In total, 157 insurance-related specialists from three Baltic countries agreed to participate in the online survey; the majority (64% of all respondents) belong to the 26–45 year-old age group, are female (78%), and from Lithuania (60%). One important reliability feature of the research sample is the variety of respondent age groups. Accordingly, a representation of five different age groups was identified, while the majority was composed of three age groups within the age range of 18–55. Nevertheless, the analysis of results shows that sociodemographic factors do not have any significant influence on digitalization, personalization, and customization evaluation or factor groups of the insurance decision-making process in digital platforms.

#### *4.2. Factor and Correlation Analysis*

#### 4.2.1. Pre-Factor Analysis

The pre-factor analysis of four indicators was conducted using the SPPS statistical software (version 26). It was oriented to investigate an internal consistency of questionnaire items, a level of questionnaire reliability, sampling adequacy, and usefulness of factor analysis. Results of the pre-factor analysis are presented in Table 4.

**Table 4.** Results of indexes in the pre-factor analysis.


As per Table 4, the Cronbach α indicator has a value of 0.875, which shows that the internal consistency of questionnaire items is relatively high and acceptable. Test reliability and acceptance to perform data reduction using EFA and CFA techniques were confirmed by the Spearman–Brown and Kaiser–Meyer–Olkin coefficients. The value of KMO was 0.839, exceeding the suggested cut-off value of 0.6. Moreover, it indicates a meritorious selection of variables. These results indicate that the sampling was adequate and supports the application of factor analysis [61]. Finally, the Bartlett's test of sphericity (a test of at least one significant correlation between two of the items studied) was significant: χ<sup>2</sup> (153) = 1140.42, *p* < 0.05. Here, the p-value is smaller than the significance level (α = 0.05); therefore, it confirms that the dataset is suitable to continue investigation within procedures of EFA and CFA.

#### 4.2.2. EFA, CFA, and Pearson's Correlation Analysis

EFA and CFA are generally accepted and widely used as a combination of statistical techniques, allowing researchers to reach high effectiveness in both statistical data content analysis and testing of construct and criteria [58]. Firstly, EFA was applied to determine a structure of latent dimensions among the observed variables reflected in the items of an instrument. Construct validity was determined by using the principal component analysis (PCA) extraction method and varimax rotation. In addition to the indicators of EFA, which were presented in the pre-factor analysis, as per Table 2, it is important to "ground" the size of the sample. Authors followed the position offered by Guadagnoli and Velicer (1988), who claim that 150 participants could be interpreted as adequate if factor loads of several items exceed 0.80 [58,62].

The CFA principal procedure was used for a better interpretation of factor structures, to test the validity of the structure obtained after EFA. Here, three groups of model fit indices were used to examine the goodness-of-fit of the model with a given dataset: videlicet the absolute fit indices (coefficients of standardized root mean square residual, SRMR), parsimonious indices (a coefficient of root means square error of approximation, RMSEA), and comparative indices (coefficients of comparative fit index, CFI; non-normed fit index, also known as the Tucker–Lewis index, TLI-NNFI) [58]. Results of these indices are provided in Table 5.


**Table 5.** Results of CFA. Source: composed by the authors, using IBM SPSS Statistics 26 (Armonk, NY: IBM Corp.).

It is important to note that the chi-square/df ratio was 1.66, indicating that the model fitted the data well. Overall, the results of fit indices in Table 5 are under the recommended cut-off points for a good model-data fit. The value of SRMR was less than 0.08, the value of RMSEA did not exceed the acceptable limit of 0.08, and values of CFI and TLI-NNFI indices were above the recommended value of 0.90. Finally, the following five groups of related factors were identified after completing EFA and CFA:


Results of the statistical factor analysis show the existence of five-factor groups, but only three of them are directly related to the digital behavior and purchase decision-making process of Baltic insurance customers: the first-factor group, the second-factor group, and the fourth-factor group. Additionally, calculation of the total rank-sum was conducted to identify the most influential factor groups; the first-factor group with 6648 points. The second-factor and the fourth-factor groups have rank-sums of 4482 and 4352 points, respectively. The Pearson correlations between the five-factor groups were calculated in a parallel way. Details of the Pearson's correlation calculation are provided in Table 6.

It should be outlined that three factors groups are identified as having a strong positive Pearson correlation: the first-factor group, the second-factor group, and the fourth-factor group. This result has a two-fold outcome. First, it supports the previous finding of these three factor groups being predominant factors within insurance purchase–decisionmaking processes at Baltics online insurance platforms. Second, it grounds the theoretical assumption about the digital environment and technological factor influences to end-users and extends the foundation of the theoretical insurance decision-making framework, with new components in the pre-purchase and purchase stages.


**Table 6.** The calculation of Pearson's correlation. Source: composed by authors by using using IBM SPSS Statistics 26 (Armonk, NY: IBM Corp.).

\*\* Correlation is significant at the 0.01 level (2-tailed). \* Correlation is significant at the 0.05 level (2-tailed).

Finally, the aim was to analyze the views of insurance-related specialists regarding their sociodemographic characteristics, i.e. age, country, and gender. As the normal distribution assumptions were not met, Mann–Whitney U and Kruskal–Wallis tests were applied to identify whether there were any statistically significant differences between different groups of respondents in the three most influential factor groups. Results of the Kruskal–Wallis test are presented in Table 7.

**Table 7.** Results of Kruskal–Wallis H test. Source: composed by the authors by using IBM SPSS Statistics 26 (Armonk, NY: IBM Corp.).



**Table 7.** *Cont.*

Results of the Mann–Whitney U test are presented in Table 8.

**Table 8.** Results of Mann–Whitney U test. Source: composed by the authors by using IBM SPSS Statistics 26 (Armonk, NY: IBM Corp.).


The Kruskal–Wallis identified that the difference between the mean values of the second-factor group for age groups was statistically significant (*p*-value = 0.004 is less than the significance level 0.05). It was identified that the highest mean rank of 103.53 was in the age group 46–55, compared to the lowest mean rank of 59.22 in the age group 18–25. Moreover, it was noticed that the value of a mean rank increases constantly within age groups. Several conclusions can be formulated in regards to the above-defined observations:


The Mann–Whitney U test was applied to test for gender differences in evaluations of factors groups. Results, presented in Table 8, do not show any statistically significant differences in the evaluation of factor groups according to gender.

#### **5. Discussion and Limitations**

First, the above-presented factor and correlation analysis validate theoretical assumptions of internal and external main driving factors and the framework of the digital insurance decision-making process, and disclose the actual Baltic insurer's perception towards this process reflection in digital distribution platforms. The main finding of a strong, positive Pearson correlation among three factor groups both support and extend the findings of Milner and Rosenstreich (2013) [30] and Rocha and Botelho (2018) [63] studies. The personal insurance experience factor has the highest rank (7.4) among all influential factors, supporting the findings by the structural model of attitude towards insurance (ATI) by Rocha and Botelho (2018) [63], where factors of perception of risk in relation to the good/asset and personal concern with finances are indicated as having the highest positive influence towards consumers' willingness to pay for insurance. Additionally, the high rank (7.2) and influence of the insurance service provider brand factor confirm marketing

domain-oriented findings of Milner and Rosenstreich (2013) [30], where marketing mixes are recognized as an outlying component of the final financial decision-making frameworks, and the structural ATI model of Rocha and Botelho (2018) [63], where the factor of trust in the industry is also outlined among three of the most influential factors. From the perspective of the Baltic insurance market research, the empirical investigation and results of marketing-oriented factors extend the finding by Kiyak and Pranckeviˇciut¯ e (2014) [ ˙ 11], where no significant statistical relationship has been identified between the insurance purchase intention and the discount factor. Moreover, this contradiction of evaluation results in insurance research on marketing-oriented factors will reveal a possible scope of the marketing domain inclusion in further studies on digital insurance decision-making. Second, the findings on the insurance literacy factor have a two-fold meaning and scientific contribution. The lower average evaluation rank (6.9) on the insurance literacy factor partially supports the findings by Weedige and Ouyang [18], Weedige et al. [64], and Allodi et al. [17], where insurance literacy was identified as a potential mitigation action and an influential factor to reduce the problem of underinsurance and low-level insurance knowledge. Otherwise, recognition of the insurance literacy factor in non-standard combination, with the factors product terms and conditions acceptability, customization, and personalization foster an additional scientific discussion and contribute, as a standpoint, to further conceptual modeling and empirical analysis.

The research also has methodological and empirical limitations. First, in methodological terms, the validity and reliability of CFA outcomes within the selected sample of 157 respondents are questionable. It contradicts the Hoelter's critical N (C.N.) recommendation, where the acceptable size of the sample (N) should be above 200 to accept a model at the 0.05 level [65,66]. Nevertheless, authors follow the methodological position by Guadagnoli and Velicer (1988) [62], claiming that a sample's size, as a function of the number of variables, is not an important factor to determine the sample pattern as being a stable and approximate population pattern. It is outlined that a research sample size between 100 and 200 observations is valid, and under this sample, "a correlation coefficient becomes an adequate estimator of the population correlation coefficient" [62]. Additionally, the value of the KMO of 0.839 confirms the adequacy of sampling. Additionally, the scope of the research instrument and assessment scale, by including more than 10 items, is considered to bring some negative effects on the reliability of measurement results and affect the perception of respondents [67]. Otherwise, it is confirmed that if the index of the Cronbach's alpha is equal to 0.7 or above, the research instrument of 10 or more items is treated as sufficient [61]. Previous research in this field also confirmed that 9-point scale scores correlated better than the 5-point scale, and, overall, an increase in the reliability level was noticed due to the transition to a higher number scale [67]. The main empirical limitation of the research can be interpreted as a future research discourse within this topic. Application of the traditional fuzzy analytic hierarchy process (FAHP), the modified version of extended fuzzy analytic hierarchy process (E-FAHP), or a combination of different multi-criteria decision-making (MCDM) techniques can reduce the vagueness inherent in a stand-alone linguistic assessment and an uncertain comparison of opinions [68]. Accordingly, future scientific discussions and empirical investigations should be focused on using this set of methodological techniques in continuous scientific investigations. A more extensive and overall continuous validation of the authors' suggested influential driving factors within a higher sample of external customers in the Baltic insurance market is required. Additionally, a simultaneous theoretical analysis within the presented topic, by adding the foundation of organization buyer behavior (OBB) and/or technology acceptance models (TAM), should be completed by the authors.

#### **6. Conclusions**

The theoretical analysis of scientific literature and recent practical trends have revealed that the global insurance market and insurers are in a multidimensional transition period, including an intensive digital acceleration (mostly due to COVID-19), dynamic changes

in the market structure, and a spread of the new hybrid operational business model. The digital environment, big data analytics, and technological and cognitive-emotional factors have become major drivers of the modern insurance concept. Otherwise, it has also been recognized that insurance customer behavior and decision-making processes in digital platforms reflect on the fundamental outcomes of the traditional CDM and hybrid choice models. Parallelly, the following features of the complex and simplified customer behavioral models, expected utility theory and behavioral reasoning theory, could be identified: interactive multi-step logic, the complexity of perceived risks, personal bias, and contextual variables in decision-making, influence of the marketing domain, and feedbackrepurchase loop. The conceptual modeling confirmed that, within digital platforms, the behavior and purchase decision making of insurance customers is a continuous, but not a simultaneous, sequence of three-stage processes. Additionally, the statistical factor analysis confirmed that the combination of three internal and external factor groups influence the insurance purchase–decision-making in Baltics digital insurance platforms. Here, a strong, positive Pearson correlation was identified among six internal factors of personal evaluation and considerations, four external factors of technological features and marketing domains, and four combined factors of insurance knowledge, levels of product customization, and service personalization. Otherwise, the analysis of sociodemographic factors, such as age group, gender, and country, showed that there are no statistically significant relationships among these types of factors and evaluations of digitalization, personalization, and customization levels within the Baltic non-life insurance market and digital insurance platforms.

**Author Contributions:** Conceptualization, G.B.; methodology, G.B. and A.G.R.; software, G.B. and A.G.R.; validation, G.B. and A.G.R.; formal analysis, G.B.; investigation, G.B.; resources, G.B.; data curation, G.B.; writing—original draft preparation, G.B.; writing—review and editing, G.B. and A.G.R.; visualization, G.B.; supervision, A.G.R. Both authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The data used in this research work are available from the corresponding authors upon request.

**Acknowledgments:** Technical software IBM SPSS Statistics 26 (Armonk, NY: IBM Corp.) support was provided by a colleague Ruta Juozaitien ¯ e.˙

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-4510-3