Bi-Objective Optimization for Multi-Row Facility Layout Problem Integrating Automated Guided Vehicle Path

For machining workshops that use automated guided vehicles (AGVs) for material handling and management, workstation layout and AGV path are coupling factors affecting the material handling cost (MHC). The multi-row layout is a typical pattern in many machining workshops, and currently, there is a lack of studies on the multi-row layout problem (MRLP) while taking into account the AGV path. This study established a bi-objective multi-row layout optimization methodology integrating the AGV path to minimize MHC and area occupancy. Specifically, workstations and transfer stations were arranged in the workshop following several non-intersecting AGV paths to decrease the material handling distance between workstations. First, a multi-row layout optimization model was established. Second, a hybrid approach combining the non-dominated sorting genetic algorithm-II (NSGA-II) and tabu search (TS) was proposed to solve it. The effectiveness of the proposed model was verified in the practice of a structural components machining workshop, and the results were compared with that of a loop-based layout. In addition, the proposed approach was compared with an exact approach and another hybrid approach based on a genetic algorithm (GA) and simulated annealing (SA). The experimental results showed that the proposed approach was able to achieve better sets of Pareto solutions within reasonable computational time.


I. INTRODUCTION
The Facility layout problem (FLP) refers to the problem of optimizing the material handling cost (MHC) or other objectives determining the most efficient arrangement of given elements on the factory floor [1]. For manufacturing factories, 20-50% of the total operating cost and 15-70% of the total manufacturing cost are attributed to MHC [2]. An optimal layout can reduce the cost of the manufacturing plant by 10-30% [3] and significantly improve productivity.
There are three common layouts used in manufacturing workshops according to the arrangement shape and the material handling path, and they are single-row, double-row, and multi-row layouts.
This study focuses on MRLP, which is often employed in machining workshops.
With the rapid development of intelligent technologies, large quantities of automatic guided vehicles (AGVs) have been used in industries for material handling to achieve automated manufacturing [21], [22] and improve efficiency and productivity. In those workshops, workstation arrangement, transfer station settings, and AGV path planning are coupling factors that need to be considered in MRLP. For example, the workstation arrangement can influence the AGV path, and the AGV path and transfer station settings can also affect the workstation arrangement.
Unfortunately, existing MRLP models did not take the AGV path into account. Optimization of facility layout and AGV path are carried out independently, which will not lead to overall optimized operation on the workshop level.
Therefore, in this study, we propose an optimization method that considers the AGV path, namely IAMRLP. In IAMRLP, the workstations and AGV transfer stations are arranged into multiple rows to ensure optimum execution of transport of materials and work pieces. The transfer stations divide the AGV path into several non-intersecting segments to avoid collision and deadlocks. Each segment has one bidirectional AGV to handle the logistics within it. The cross-segment logistics are first disposed to the transfer station and then transported to the destination using the AGV within the segment. The general problem is known to be NP-hard [23]. To solve the problem in an acceptable time, a hybrid approach combining the non-dominated sorting genetic algorithm-II (NSGA-II) with tabu search (TS) is adopted. NSAG-II is used to effectively determine the sequence and location of each workstation, and TS is used to optimize the Pareto solutions obtained by NSGA-II to avoid falling into a local optimum. The proposed method is verified through a case study and numerical experiment.
Key contributions of this paper are as follows: 1) It proposes a new machine layout problem considering the AGV path called an IAMRLP.
2) It establishes a Bi-objective optimization model for this new problem.
3) It applies the proposed approach to a structural machining workshop.
The remainder of this paper is organized as follows: Section II describes the background and research status of MRLPs. Section III presents the optimization problem in this study and the proposed model formulation. The hybrid solution approach is described in Section IV. Section V presents a case study of the proposed method, and a numerical experiment is performed in Section VI. Lastly, the conclusions are summarized in Section VII.

II. LITERATURE REVIEW
MRLPs are reviewed below since they are closely related to our work.
MRLP is a class of facility layout problems that determine the arrangement of facilities in certain rows to minimize MHC [24]. Research on MRLP was taken place since the early 1960s [25]. Koopmans and Beckmann modeled MRLP as a quadratic assignment problem [26]. Kusiak and Heragu surveyed models and algorithms for MRLP [27]. Heragu and Kusiak proposed a triangle assignment algorithm to solve multi-row cluster machine layout in a flexible manufacturing system (FMS) [28]. Fischer et al. [29] presented new mixed-integer linear programming models for the MRLP. Sooncharoen et al. [30] presented a novel biogeography-based optimization tool to solve the unequal area MRLP to minimize the total material flow distance. Ahmadi et al. [31] carried out a literature review on MRLP, and presented the applications, essential features, approaches, and resolution methodologies on MRLP.
Liu et al. [32] established a multi-row mix integer programming model for FMS, aiming to minimize the total cost of logistics handling. A GA was used to analyze the problem. Forghani et al. [18] presented a new approach for multi-row cellular manufacturing system layout design to minimize the total material handling cost. A four-stage heuristic method was proposed to solve the model. Safarzadeh and Koosha [24] proposed an extended MRLP model with fuzzy constraints. The objective function was to minimize the material handling and lost opportunity costs, and a GA was used to find the best solutions. Tubaileh and Siam [16] addressed the problems of locating machines in single, double, and multi-row layouts for FMS to optimize the total cost of material transportation. The best machine arrangements were obtained using ant colony and simulated annealing algorithms. Hu and Yang [33] proposed a mathematical MRLP model in semiconductor fabrication plants by minimizing the total transportation distance, and a particle swarm optimization algorithm was used to solve the model. Wan et al. [20] established an MRLP model with extra clearances, with the objectives of minimizing material flow cost and layout area. A hybrid approach combining an improved multi-objective greedy randomized adaptive search procedure and linear programming was proposed for the problem. Anjos and Vieira [1] proposed a new mixed integer linear programming model for MRLP, and a two-stage optimization algorithm was used for large instances solutions. Herran et al. [19] proposed a variable neighborhood search algorithm to solve the multi-row facility layout problem.
The AGVs are increasingly used in modern workshops. The facility layout and AGV path planning are coupling factors for those workshops. More AGVs are required, and transport control is more complicated in case of a bad layout.  Unfortunately, all the above MRLPs only considered the arrangement of workstations, and the AGV paths were not considered. Therefore, they are difficult to reach the optimum states. To the best of our knowledge, there is no report on MRLP considering AGV paths.

III. MODEL FORMULATION A. PROBLEM DESCRIPTION
Among the common layout designs, the multi-row layout has the advantages of high area occupancy, high flexibility of logistics paths, and short logistics distribution distances. Hence, it has been widely used in manufacturing workshops. Therefore, this study uses the multi-row layout problem as the research subject.
For the multi-row layout, the materials are distributed by AGVs along lateral and horizontal paths. There are four common types of AGV systems: network system, tandem system, single loop system, and segmented system, as shown in Fig. 1. The network AGV system ( Fig. 1A) has good flexibility, yet the AGV paths often overlap, resulting in collisions, blockages, and deadlocks. The single loop system sets the AGV path along the shortest route, and the AGVs run unidirectionally (Fig. 1B). The segmented system sets the AGV path along the shortest route. However, the paths are divided into several non-intersecting segments that are connected by transfer stations (Fig. 1C). The AGVs run bidirectionally within each segment, the cross-segment logistics is first distributed to the transfer station. And then transported to the destination using its inner-segment AGVs. The tandem system divides the workshop into several non-interleaved areas, transfer stations are set in each area to connect with other areas, the AGV only operates within the area, and cross-area logistics are completed through transfer stations (Fig. 1D).
Single loop, segmented, and tandem systems can avoid AGV collisions, blockages, and deadlocks. However, they still have some limitations. For instance, for single loop and segmented systems, limited by the path and direction, the AGV travel distance per task is longer than that of other systems. In addition, for a tandem system, the area occupancy is larger than other systems due to the use of transfer stations.
To simultaneously reduce the length of the logistics path and area occupancy, an optimization model for MRLP integrating the AGV path is proposed, as shown in Fig. 2. In the model, the AGV path is divided into several segments by the transfer stations, and there is a route from each transfer station to the path to which it belongs, aiming to reduce the logistics distance. In each segment, an AGV is running bi-directionally. When the AGV needs to distribute material between segments, the materials are first distributed to the transfer station and then transported to the destination. For example, to distribute materials from workstation 8 to workstation 3, the AGV in segment two first distributes the materials to transfer station T1 between segments 1 and 2. Then, the AGV in segment 1 distributes the materials from T1 to workstation 3.

B. ASSUMPTIONS
The following assumptions are made for the MRLP model incorporating the AGV path.
1) The workstations are rectangles with known lengths and widths.
2) AGVs only move horizontally or laterally in the workshop.
3) The feeding and blanking position of each workstation is the same.
4) The workstation is placed in parallel with the workshop. 5) When the x-axis of a workstation exceeds the workshop's width, it will be placed on the next row.
6) The transfer stations can only be set between adjacent workstations, and their sizes are ignored. 7) When workstations and transfer stations are in the same row, their ordinates are the same.

C. PARAMETER VARIABLES
In Fig. 2, the x-axis is the horizontal direction, the y-axis is the lateral direction, L is the length of the workshop, W is the width of the workshop, T i is transfer station i, l i is the length of workstation i, and w i is the width of workstation i.
x i is the abscissa of workstation i, and it is calculated as follows.
When workstation i is in odd rows: When workstation i is in even rows: Especially, when i is the first workstation in an odd row: When i is the last workstation in an even row: y i is the ordinate of workstation i, and it is calculated as follows: For workstations in the first row: For workstations on other rows: where the set Rn is the set of workstations in the n th row, W is the width of the workshop. dx i is the minimum safety distance of workstation i in the horizontal direction. dy i is the minimum safety distance of workstation i in the lateral direction.
x tp is the abscissa of the p th transfer station. When the p th transfer station is behind workstation i in the same row, x tp is calculated as follows: y tp is the ordinate of the p th transfer station. When the p th transfer station and workstation i are in the same row, y tp is calculated as follows: dt i1 is the distance between workstation i and the first transfer station on its logistics path, dt i1 is calculated as follows: dt pj is the distance between the last transfer station and workstation j on its logistics path: dt i t i+1 is the distance between transfer station i and transfer station(i + 1): dw i w i+1 is the distance between adjacent workstations: d ij is the distance between workstations i and j, unit: m. d ij is calculated as follows: When there is no transfer station between workstations i and j: When there are transfer stations between workstations i and j: Here, P denotes the number of transfer stations.

D. CONSTRAINTS
The MRLP model in this study is constrained by the workstation size, workshop size, etc. The constraints are as follows: where N is the number of workstations, R is the number of rows, G ir = 1 means that workstation i is in row r, G ir = 0 indicates that workstation i is not in row r.
Equations (15) and (17) mean that the layout cannot exceed the size of the workshop, the horizontal spacing of each row is no greater than the length of the workshop, and the lateral spacing of all rows is no greater than the width of the workshop.
2) SPACING CONSTRAINTS adjacent workstations in the same row do not overlap, the minimum horizontal safety distance must be satisfied, as expressed in (18), and the spacing of workstations between adjacent rows must satisfy the minimum lateral safety distance, as expressed in (19).
3) EQUATION (20) (21) where, F is the total MHC between the workstations, unit: kg · m; N is the total number of workstations to be arranged, q ij is the material flow from workstation i to j, unit: kg. Objective Function 2: Minimizing the area occupancy: where, S is the area occupancy, unit: m 2 .
Hence, it has become a popular algorithm for solving FLPs. In the early iterations of NSGA-II, the search speed is fast. However, as the number of iterations increases, the solutions are gradually homogenized because they mainly depend on mutation. TS is a global optimization algorithm that can obtain excellent neighborhood solutions. Hence, it is expected to improve the overall performance by integrating NSGA-II with TS. Therefore, a hybrid algorithm integrating NSGA-II and tabu search, namely NSGA-TS, is proposed in this study. NSGA-II is used to obtain the initial solution efficiently, and TS is applied to obtain reasonable neighborhood solutions to prevent falling into a local optimum. The flow chart of NSGA-TS is shown in Fig. 3.
Step 1: The workstations are coded using real numbers, and the transfer stations are coded using binary numbers.
Step 2: Generate the initial population, set the population number POP and the number of iterations GEN , set gen = 1, and set the tabu list empty.
Step 3: Calculate the objective value of each individual, perform rapid non-dominant sorting, and calculate the congestion distance of the contemporary population.
Step 4: Individuals are selected by linear sorting and then crossed and mutated to generate the child population.  Step 5: Combine the parent and child populations, use the elite retention strategy to generate the new population, perform rapid non-dominant sorting, and calculate the congestion distance of the new population.
Step 7: Use the optimal TS solution to replace the most dominant individual in the new population, set gen = gen+1, and continue until the number GEN is reached.

A. ENCODING AND DECODING
In this study, the arrangement of workstations and transfer stations needs to be defined first. We use a two-stage method to encode the chromosomes. The coding method is illustrated in Fig. 4.
The first half of the chromosome represents the order of the workstations and is coded using natural numbers that represent the workstation number in set N . The second half of the chromosome represents whether there is a transfer station behind the workstation corresponding to the exact location of the first half, which uses binary encoding. ''1'' indicates that there is a transfer station, and ''0'' indicates that there is not a transfer station. According to Fig. 4, the workstation layout is [1,5,4,2,7,3,8,9,6], and the transfer station layout is [0, 0, 0, 0, 1, 0, 1, 0, 0], that is, the workstations are placed according to the gene chain, and there is a transfer station behind the workstations No. 7 and No. 8. The decoding process is as follows: Step 1: Arrange the workstations by order of the chromosome and their coordinates.
Step 2: Arrange the transfer stations according to the second half of the chromosome.
Step 3: Manually connect each transfer station to the path to section the AGV paths into the segments.

B. POPULATION INITIALIZATION
In this study, the population initialization is implemented in three steps: First, generate the initial population randomly. Calculate the coordinates of each workstation and transfer station and the distance between workstations.
Then, calculate the values of the objective functions. Finally, perform rapid non-dominant sorting and calculate the congestion distance of the contemporary population.

C. SELECTION
To retain individuals with high fitness values in the next generation so as to improve the convergence and efficiency of the algorithm, after non-dominated sorting and crowding calculation of the population, the linear sorting method is used to select the individuals that need to cross and mutate in the parent generation. The individuals are first sorted according to their non-dominated sorting value. A lower sorting value is preferable to a higher value. A higher congestion distance is preferable if the sorting value is the same.
In the linear sorting method, all individuals are numbered. The best individual number is n, the probability of being selected is p max , the worst individual number is 1, the probability of being selected is p min , and the probability of selecting other individuals is based on the following equation: In this study, p max = 0.9, p min = 0.1.

D. CROSSOVER
The crossover operation is performed on different chromosome sections to avoid invalid solutions. The distance between workstations is the main factor when calculating MHC. The smaller the distance, the smaller the MHC. Therefore, we propose a heuristic crossover operator based on the nearest distance between workstations. In the crossover operation, if the current gene of the offspring chromosome is O ′ i1 , then the subsequent workstations of O ′ i1 in the parent chromosome are P ′ i1 and P ′ i2 . The workstation closer to O ′ i1 is selected as the subsequent gene O ′ i2 of the offspring chromosome. This operation is repeated until the end of the parent chromosome. The heuristic crossover strategy can extract genes in the parent generation to optimize the offspring chromosomes. The strategy is implemented in 8 steps: Step 1: Select two parent chromosomes, P 1 and P 2 , to perform the crossover operation.
Step 2: Select the first gene of P 1 and put it into chromosome O ′ ,and delete the same gene from P 1 and P 2 , as shown in Fig. 5.
Step 3: Take the last gene in O ′ as O ′ i1 , find the genes right to O ′ i1 from P 1 and P 2 , as P ′ i1 and P ′ i2 , calculate the distances between workstation O ′ i1 , P ′ 1 and O ′ i1 , P ′ 2 , as dist1 and dist2, as shown in Fig. 6.    Step 4: If dist1 < dist2, place P ′ 1 into chromosome O ′ , delete P ′ 1 from P 1 and P 2 ; otherwise place P ′ 2 into chromosome O ′ , and delete P ′ 2 from P 1 and P 2 , as shown in Fig. 7.
Step 5: Repeat Step2 to Step4 until half the length of O ′ i1 is filled, as shown in Fig. 8.
Step 6: Combine chromosome O ′ with the remaining genes of P 1 and P 2 to generate individual O 1 and O 2 , as shown in Fig. 9.
Step 7: Repeat steps 1-6 on P 2 to generate individual O 3 and O 4 .
Step 8: Select the two best individuals in O 1 , O 2 , O 3 , and O 4 as the offspring chromosomes.
A two-point crossover method is used for the transfer station section. First, the start and end points are randomly selected, and the chromosomes between the start and end points are exchanged. Conflict detection is performed  according to assumption No. 6, where the transfer station can only be between adjacent workstations. That is, in the transfer station section, a gene with code 1 cannot be next to another gene with code 1; if so, the second gene 1 is changed to 0.

E. MUTATION
Mutation operation uses the neighborhood variation method. First, generate an array of random binary numbers (0 and 1) with a length 2n, where n is the length of the chromosome. Second, the chromosome is mutated at the first element ''1'' in the array. For workstations, the corresponding genes are exchanged. For transfer stations, the gene is changed to the opposite value, and conflict detection is performed according to assumption No. 6. The mutation process is shown in Fig. 10.

F. TABU SEARCH
To prevent a premature algorithm, TS is added to the algorithm. TS is performed using the following steps: Step 1: Take the individual with the lowest dominance level as the initial solution of TS, and initialize the tabu list.
Step 2: Perform a neighborhood search to obtain candidate solutions. In this study, we use the insert operator to generate neighborhood solutions. Neighborhood solutions are generated by removing one workstation from the chromosome and inserting it at a different position on the same chromosome. First, randomly select a number p from 1 to m, m is half of the number of workstations n. Then, the p th gene of the workstation section is removed and inserted at a different position. The insertion operator is only performed for the workstation section while the transfer station section remains unchanged. For example, suppose the randomly chosen number is 3. In that case, the corresponding gene is first removed  from the chromosome, as shown in Fig. 11, and then inserted behind each of the remaining genes to generate neighborhood solutions. Only the workstation section performs the insertion operator, and the transfer station section remains unchanged.
Step 3: Perform the rapid non-dominant sorting and calculate the congestion distance.
Step 4: When the candidate solution with the lowest dominance is better than the initial optimal solution, replace the initial solution with the candidate solution. Otherwise, add the candidate solution to the tabu list. Then, repeat TS until the termination condition is met.

V. CASE STUDY
In order to verify the MRLP model proposed in this study, a structural components machining workshop in Hubei, China, is used as a case study. The workshop belongs to the customer-oriented discrete manufacturing mode. Since the workshop customers are relatively fixed, the product quantity of each plan period mostly stays the same. The logistics quantity in Table 1 is the daily average over six planning periods.
The length of the workshop is 42m, and the width is 30m. The workshop has 22 machines and 16 production types. Table 2 shows the size of each machine. According to the actual situation, the minimum lateral safety distance between adjacent workstations is 2m, and the minimum longitudinal safety distance between rows is 2m.
The original facility layout is based on the cluster principle, as shown in Fig. 12. In this layout, there is a high MHC and low production efficiency. Therefore, our proposed approach is used to optimize the layout. The parameters are as follows: the initial population size of the NSGA-II algorithm is 200, the number of generations is 200, the maximum probability of being selected p max is 0.9, the minimum probability of being selected p min is 0.1, and the TS iterations is 50.
Pareto solutions obtained by our proposed approach are shown in Fig. 13. The horizontal axis represents the MHC, the  vertical axis represents the area occupancy, the 'o' represents the original layout, and the ' * ' represents the Pareto solutions. As shown in Fig. 13 and Table 3, the NSGA-TS obtained better solutions than the original layout. The MHC and the area  occupancy were reduced by 10.1% and 4.2%, respectively. When we increased the number of generations to 10000 to test whether it could improve convergence, we observed that the Pareto solutions did not improve with the number of generations. Therefore, we set the number of generations in the numerical experiment to 200.

VI. NUMERICAL EXPERIMENT
In this section, we perform a numerical experiment for the proposed model and the NSGA-TS. Since IAMRLP is a new problem proposed in this paper, there exist no solution approaches for it, and hence, our approach is compared with an exact approach (CPLEX) and GA-SA, which is an effective and popular method used to solve FLP [39], [40], [41]. To verify the effectiveness of our approach, the result is compared with a loop-based facility layout.

A. PROBLEM SET AND PARAMETER SETTINGS
The workshop data in Section V are used as the problem set P1. To test the ability of NSGA-TS to handle large-scale  problem instances, we duplicate the workstation size data ten times to build a large-scale problem set P2 that consists of 220 workstations. The logistics quantity is randomly generated, and the workshop size is expanded to 150m ×90m. The procedure proposed in Forghani et al. [40] and the encoding, selection, crossover, and mutation operator proposed in Section IV are used to build the GA-SA. The other parameters are listed in Table 4.
To solve the formulation in Section III, the two objectives (21) and (22) are combined linearly to form a single objective as shown in (26).
In this study, α increases from 0 to 1 by a small step size of 0.05 to find more Pareto solutions.

B. RESULTS AND COMPARISONS
Five independent runs are performed for each problem instance.

1) PROBLEM SET P1
Pareto solutions are shown in Fig. 14. Table 5 shows the max, min, and mean values of MHC and area occupancy, and the mean runtime over the five runs. Solutions obtained by NSGA-TS are better than those of GA-SA. Since the FLP problem has a high computational cost, the exact approach CPLEX only obtains some Pareto solutions within the runtime, and all are inferior to those of NSGA-TS. The runtime of NSGA-TS and GA-SA is close.

2) LARGE-SCALE PROBLEM SET P2
Pareto solutions are shown in Fig. 15. Table 5 shows the max, min, and mean values of MHC and area occupancy, and the mean runtime over the five runs.
The result shows that for a problem with a size of up to 220, although NSGA-TS takes a longer runtime than GA-SA, the Pareto solutions of NSGA-TS are better than those of GA-SA, IAMRLP decreased the MHC by 12%. Within the given runtime, CPLEX does not obtain any Pareto solution at all. When the runtime is increased to 24 hours, CPLEX still cannot obtain any Pareto solution.

C. COMPARISON WITH A LOOP-BASED LAYOUT
Previous studies showed that the loop construction method provided an optimal solution in the integrated solution of FLP and AGV path planning. In this section, we compare our approach with a loop-based facility layout shown in Fig. 16.
Each inner loop has a transfer station for cross-loop transportation. An AGV is used for logistics handling. Crossloop logistics are transported between transfer stations along the outer-loop route using AGVs. The positions of workstations, the distance between workstations, and the positions of transfer stations are calculated as illustrated in [42]. The   objective functions are the same as those in (21) and (22). The algorithm used is NSGA-TS. The chromosome encoding method proposed in Section IV is used. One difference is that number 1 indicates that the corresponding position of the workstation section is the end of the loop. Table 6 shows the MHC and area occupancy's max, min, and mean values.
The Pareto solutions are presented in Fig. 17. The results show that IAMRLP leads to a lower MHC and a smaller area occupancy. IAMRLP reduced the MHC by 4% and the area occupancy by 6.3%.

VII. CONCLUSION
In this study, an IAMPLP is proposed. Different from traditional MRLPs, this problem takes into account AGV paths. The main contribution of this study is integrating AGV paths into the MRLP model to optimize material transportation and AGV routing between workstations. A bi-objective model is established to identify the workstation's and the transfer station's sequence on multiple rows and their exact locations. A hybrid algorithm is used to solve it. The proposed IAMPLP can be used to optimize facility layout and AGV path planning and help manufacturing workshops achieve high efficiency and productivity.
In the next step, we will consider the influence of changing material flow between different planning periods and the rearrangement costs, i.e., a dynamic facility layout problem (DFLP), and we will integrate digital-twin with heuristic algorithms to solve the FLP problem while considering the actual equipment geometry, so as to achieve a more realistic facility layout.