Heuristic Approach for a Combined Transfer Line Balancing and Buffer Allocation Problem Considering Uncertain Demand

Featured Application: This research was initiated by an industrial project. The problem was the design and conﬁguration of machining lines for engine blocks. The proposed approach was validated using four real cases provided by the industrial partners of the project. The proposed approach could easily be applied to the design and conﬁguration of any machining line for the production of a single complex mechanical component. Abstract: In this paper, we refer to a real case study of an industrial partner recently committed to its project on the design of a multi-unit and multi-product manufacturing system. Although the considered problem refers to an actual complex manufacturing system, it can be theoretically classiﬁed as a union of two key problems that need to be solved during the transfer line design stage: the transfer line balancing problem (TLBP) and the buffer allocation problem (BAP). As two closely related problems, TLBP and BAP usually have similar optimizing directions and share the same purpose: ﬁnding a balance between the performance of the transfer line system as well as minimizing investment costs. These problems are usually solved sequentially, but this leads to solutions close to a local optimum in the solution space and not to the global optimum of the overall problem. This paper presents a multi-objective optimization for concurrently solving transfer line balancing and buffer allocation problems. The new approach is based on a combination of evolutionary and heuristic-based algorithms and takes into account the uncertainty of market demand. To validate the proposed approach, an industrial case study in a multi-unit manufacturing system producing multiple products (four engine blocks) is discussed.


Introduction
Transfer lines are serial machining systems dedicated to the production of large series, which are widely used production systems in machining environments. They are composed of a set of workstations and an automatic handling system. Each workstation carries out one identical set of operations every cycle. The design of transfer lines is comprised of several steps: product analysis, process planning, line configuration, transport system design, and line implementation [1]. During the whole design stage, there are two important research issues: the transfer line balancing problem (TLBP) and the buffer allocation problem (BAP), which belong to the two different steps, "line configuration" and "transport system design", respectively.
As two closely related problems, TLBP and BAP usually have the similar optimizing directions and share the same purpose: finding a balance among the system performance of the transfer line in relation to minimizing investment costs, maximizing production efficiency, etc. With the further study of the TLBP and BAP, researchers from both sides have generally shown a growing interest in simultaneously solving the two problems considering their closeness.
From the line balancing side, Battaïa et al. [2] presented a more recent and relevant survey, where they focused on the articles solving these aggregated problems that involved both the balancing problem and its "neighbor" problems during the design stage of the production line. The authors classified and discussed three examples of problem combinations, which were "the process selection and line balancing", "the line layout design and balancing", and "the line balancing and task sequencing". Finally, it was concluded that "it is necessary to consider several decision problems simultaneously to enhance the final line performance".
From the buffer allocation side, in a recent review on BAP in production lines [3], Weiss et al. concluded that "The buffer allocation is only one among several design decisions, which include the decision on the number of machines per station, the characteristics of the machines per station, and the workload allocation among the stations. In the literature, these decisions are typically optimized individually. However, in the design phase of a real flow line, a joint decision is often possible." Therefore, solving the two decision problems simultaneously in order to enhance the final performance of a machining line, and to thereby become closer to the real industrial environment seems to be a promising path.
In this paper, we refer to a real case study of an industrial partner recently committed to a project on the design of a multi-unit and multi-product manufacturing system. Although the considered problem refers to an actual complex manufacturing system, it can be theoretically classified as a union of the two key problems TLBP and BAP. The factory is divided into four manufacturing units, where each unit is in charge of the machining of one of four different engine blocks. Each unit may consist of a number of stations with parallel machine tools and different station configurations (part orientation, fixture, and datum system). Moreover, there is a buffer area between two stations that decides the buffer capacity. From the perspective of the nature of the problem, in order to design this kind of manufacturing system, all the detailed decisions inside the manufacturing units have to be made, which includes the selection of the number of stations, the fixturing type per station, the number and type of machine tools per station, the operations allocated to each station, and the buffer allocation between two stations. Meanwhile, after this series of sub-problems are solved, the selected sub-solutions will work together and become coupled to the overall system on different levels, thus greatly affecting the final performance of the manufacturing system. As the focus of decision makers, the final system performance can be quantified into a set of evaluation criteria, such as production efficiency, investment cost, floor space, etc. Thus, when solving this kind of complex problem, the industrial company has shown great interest in choosing between a set of feasible solutions with different final performances, instead of solving sub-problems by breaking down the complete problem.
Therefore, we propose to solve TLBP and BAP simultaneously with a multi-objective heuristic approach. We also refer to the real manufacturing system in the automotive sector with four units producing four different kinds of engine blocks to further validate the proposed approach and compare the results on all the four real situations. The four kinds of engine blocks are shown in Figure 1, where part A is the same case study presented in [4]. The four parts have six faces and 351, 242, 421, and 285 machining tasks, respectively. They have 49, 50, 57, and 58 different machining features with 84, 68, 109, and 96 different machining operations, respectively.
The rest of this paper is structured as follows: Section 2 presents a literature review on some references from the research area of TLBP, BAP, and the simultaneous approaches to them. Then, the main contribution of this paper is discussed. Section 3 presents the detailed mathematical description of the problem. Section 4 provides the proposed approach, describing in detail the heuristic principles and tuning of the parameters of the proposed approach. Section 5 presents the industrial case study of four engine blocks to compare and validate the proposed approach. At the end, the results and outcomes are summarized in the conclusion. A detailed description of the industrial case study with four manufacturing units is reported to ensure the reproducibility of the proposed approach and the possible usage of the case study as a testbed for other possible approaches.  The rest of this paper is structured as follows: Section 2 presents a literature review on some references from the research area of TLBP, BAP, and the simultaneous approaches to them. Then, the main contribution of this paper is discussed. Section 3 presents the detailed mathematical description of the problem. Section 4 provides the proposed approach, describing in detail the heuristic principles and tuning of the parameters of the proposed approach. Section 5 presents the industrial case study of four engine blocks to compare and validate the proposed approach. At the end, the results and outcomes are summarized in the conclusion. A detailed description of the industrial case study with four manufacturing units is reported to ensure the reproducibility of the proposed approach and the possible usage of the case study as a testbed for other possible approaches.

Literature Review
Although either TLBP or BAP has been the focus of extensive publications, this paper mainly focuses on the simultaneous approach to TLBP and BAP. Therefore, this section first lists some of the relevant literature in the research field of TLBP and BAP, then reviews in detail the existing approaches to the simultaneous optimization of TLBP and BAP, and finally discusses the novel contributions of this paper.

Transfer Line Balancing Problem
Transfer lines are production systems that consist of serial stations where operations are continuously performed [1]. TLBP consists of assigning the set of indivisible units of work elements to a set of sequential workstations for certain environments (dedicated, flexible, or reconfigurable transfer line). Transfer line balancing problems were first reported by Szadkowski [5] and first presented by Dolgui et al. [6] for the dedicated transfer line. Another type of balancing problem, closer to the one considered in this paper, was first introduced by Essafi et al. [7]. They considered the flexible or reconfigurable transfer lines. The transfer line has multiple stations that consist of one or more identical machine tools operating a given set of sequenced tasks, and each station is usually equipped with single spindle machine centers. Different from classic TLBP, this has two more relevant problems inside, which are called the line configuration problem and the operation allocation problem. The line configuration

Literature Review
Although either TLBP or BAP has been the focus of extensive publications, this paper mainly focuses on the simultaneous approach to TLBP and BAP. Therefore, this section first lists some of the relevant literature in the research field of TLBP and BAP, then reviews in detail the existing approaches to the simultaneous optimization of TLBP and BAP, and finally discusses the novel contributions of this paper.

Transfer Line Balancing Problem
Transfer lines are production systems that consist of serial stations where operations are continuously performed [1]. TLBP consists of assigning the set of indivisible units of work elements to a set of sequential workstations for certain environments (dedicated, flexible, or reconfigurable transfer line). Transfer line balancing problems were first reported by Szadkowski [5] and first presented by Dolgui et al. [6] for the dedicated transfer line. Another type of balancing problem, closer to the one considered in this paper, was first introduced by Essafi et al. [7]. They considered the flexible or reconfigurable transfer lines. The transfer line has multiple stations that consist of one or more identical machine tools operating a given set of sequenced tasks, and each station is usually equipped with single spindle machine centers. Different from classic TLBP, this has two more relevant problems inside, which are called the line configuration problem and the operation allocation problem. The line configuration problem defines the number of stations and the number of parallel machine tools per station, while the allocation problem defines the operations allocated to each station and their positions in the station operating sequence.
To our knowledge, this kind of TLBP was first studied by Essafi et al. [7]. The authors introduced a mixed integer programming method for which minimizing the number of machine tools was considered as the single objective function. A greedy heuristicbased approach was then presented to solve this problem [8] with the same objective function. Later, an ant colony algorithm [9] and a greedy randomized adaptive search [10] were proposed with a single objective combining the performance of line balancing and the number of stations and of machines. A set partitioning model-based approach was proposed by Borisovsky et al. [11,12] to solve the problem. The aim was to minimize the number of machine tools. TLBP was transformed into a linear integer program based on the set partitioning model, where the operation sequencing for each station was optimized. A genetic-based approach [13] was also proposed, where the decoding heuristic was designed. This approach considered a single objective function (the total number of machine tools) introducing a value punishing the infeasible solutions. In addition, the hybrid techniques of Benders decomposition and an ant colony algorithm were proposed by Osman and Baki [14] to obtain the optimal configuration of the production line by minimizing the non-productive time.
A simultaneous approach was proposed by Delorme et al. [15] combining the "balancing problem", the "scheduling problem", and the "equipment problem", and using multi-objective optimization. In this approach, the balancing problem was used to decide the allocation of machining tasks to stations and the number of stations, the scheduling problem focused on sorting machining tasks inside stations, and the equipment problem was used to select the number of parallel machines for each station. Three sub-methods were proposed and combined to sequentially solve the corresponding sub-problems. During the optimization, the authors created a solution pool to carry the feasible plans, where new feasible plans with a lower cycle time would be obtained through combining two solutions as parallel production lines. Finally, Pareto optimizing was introduced to analyze the solution pool.

Buffer Allocation Problem
Buffer allocation problems focus on finding optimal buffer capacities to be allocated to buffer areas in the production system for a specific objective. BAP is an NP-hard combinatorial optimization problem [16], and it was first presented by Koenigsberg [17]. Demir et al. [18] provided a systematic review of the studies published since 1998. They classified the literature based on four aspects: the line topology, the objective functions, the evaluative method, and the generative method. The authors also concluded that "handling BAP in a multi-objective manner still seems an important research issue". Recently, Dolgui et al. [19] and Su et al. [20] suggested studying multi-criterion optimization problems instead of optimizing a single-objective BAP.
Since the ultimate goal of BAP is to optimize the performance of the production line, an interesting research area called the bowl phenomenon has emerged under an unpaced production environment. It was introduced by Hiller and Boling [21], while a review of the phenomenon can be found in [22]. This research area considers a relevant problem that is the simultaneous allocation of the workload and the buffer capacities in an unpaced production line. Hillier and Hiller [23] tried to approach this kind of problem for a reliable production line with no parallel stations. During their further research [24], exact numerical results were shown for small lines (no more than four stations) with exponential or Erlang distributions of processing time, while some possible heuristic methods and simulation studies were provided for other situations where the production lines had more stations or different processing time distribution. Ng et al. [25] proposed a simulation-based multi-objective optimization (MOO) approach. Four aspects of production systems were considered when solving the problem, which were throughput, work in process, average idle times, and average buffer level. It was also concluded that "the ultimate goal of this kind of studies is to propose a Pareto-based MOO methodology for studying the simultaneous effects of workload balancing and buffer allocation to real production lines with complex routings and various sources of variations." It is worth noting that most of these approaches assumed that the total amount of the workload remains certain throughout the line and can be arbitrarily divided at each station. This represents the complexity of solving TLBP, where the workload of a station is determined by a set of indivisible machining operations, and each operation is constrained more or less by different types of constraints (such as precedence constraints) during the allocation.

Simultanous Approach to TLBP and BAP
Despite the extensive literature on both the TLBP and the BAP, the simultaneous approach to the combined problems has not been fully addressed. The most relevant approaches are as follows.
Tiacci [26] studied the simultaneous optimization of the assembly line balancing problem (ALBP) and BAP in the mixed model assembly line. This problem is characterized by a stochastic task time, parallel stations, and a buffer between the stations. A genetic approach was introduced, and a parametric event-oriented simulator was built to evaluate the solutions considering a single objective function "normalized design cost". It contained the investment cost together with a penalty function punishing the infeasible plans.
Koren et al. [27] combined the BAP with the scalability problem introduced in [28], which made the full problem able to simultaneously reconfigure and rebalance the line and reallocate the capacities of the buffer area. A genetic algorithm was proposed to optimize the single objective function that is a combination of the number of machine tools and the throughput.
To our knowledge, the first attempt to address TLBP and BAP simultaneously was reported in Liu et al. [29]. A GA-PSO approach was proposed with a single objective function where the production rate was combined with the total investment cost. An industrial case study was introduced to validate the method and to compare it with that of the sequential optimization. Referring to this industrial case study, although sequential optimization may obtain a better value of the objective function of balancing problem, simultaneous optimization leads to a better overall solution. Following this result, the simultaneous optimization of TLBP and BAP with a multi-objective optimization approach was carried out in [4]. The overall problem was thus to make simultaneous decisions on all the following:

•
Line configuration: the number of stations and the number of identical machine tools for each station; • Station configuration: the selection of fixturing type and machine tool type for each station; • Operation allocation: the allocation of each operation to a station and its positioning in the machining sequence; • Buffer allocation: the allocation of the buffer capacities between two stations; • In order to solve the whole problem, several heuristic principles and evolutionary algorithms were built. A priority-based coding system was designed to make the approach suitable for the application of different kinds of algorithms. Four heuristic principles were proposed to ensure the feasibility of the solutions. In order to obtain the optimal set of solutions, Non-dominated Sorting Genetic Algorithm-II (NSGA-II) and Multi-Objective Particle Swarm Optimization (MOPSO) were introduced and fine-tuned, where the total investment cost and the throughput from simulation-based evaluation were considered as the two objective functions. Finally, an industrial case study was introduced to demonstrate the validity of the proposed approach, where the comparison results showed the benefit of NSGA-II for solving the problem.

Contribution of the Paper
The main contribution of the paper is twofold: to deal with the uncertainties of the company's forecasts of the market demand and to reduce the computational time while improving the search in the solution space. The most relevant works in this research area rely on a correct estimation of the market demand, which means the optimization process always revolves around a given value. However, the market demand in the real situation is of course uncertain, and it is more reasonable to give the decision maker a set of solutions that span a range of market demands rather than referring to a single value. To consider a range of market demands, the problem must be modified with the consequence that the solution space can expand considerably, thus increasing the difficulty in solving the overall problem.
The second issue refers to the searching speed and ability of the algorithm. Due to the inner operators from the evolution algorithms, newly generated solutions have the chance to be feasible but are obviously far from the optimum. However, all solutions of each generation, even the worst, need the simulation-based evaluation of the throughput. This could lead to a waste of computing time. In addition, it may result in the request for applying again and again the simulation-based evaluation of the throughput on the identical solutions, which actually leads to overestimations: the production rate may be improved just by running the simulation on the same solution due to the natural variability of this result. This may lead to issues related to both the local optimization of the objective function and a waste in computational time.
Therefore, we propose a combination of an evolutionary algorithm and a heuristicbased algorithm to solve the problem considering the market demand interval. Firstly, new heuristic principles in the decoding stages were designed to fulfil the industrial requirement of the expected demand interval, which not only ensured each individual could be decoded into a feasible solution, but also improved the quality of the solutions during the decoding stage, thereby improving the final performance of the algorithm. At the same time, the heuristic principles gave each solution classic evaluating criteria from the line balancing side, which also provided a direct and effective basis for swiftly screening the worse solutions among the entire population. Based on this, a pre-screening strategy is proposed to eliminate the obviously worse or previously evaluated solutions in each generation. With the help of this strategy, the algorithm may reduce invalid searches and waste of computing time and thus improve the optimizing performance. In addition, a simulation strategy is introduced and was tested to enhance stability when estimating the throughput of a single solution.

Assumptions
The assumptions of the proposed problem were as follows:

1.
A flexible transfer line produces a single part and consists of at least one station in the series and buffer areas between two consecutive stations; 2.
Multiple machining features are assigned to a single part and each of them has at least one machining operation which needs to be operated in sequence; 3.
The operation time of each machining operation is known and there may be precedencies among the machining operations; 4.
A finite number of buffers can be allocated at the buffer area; 5.
Each station has one or more machine tools that are identical and use an identical type of fixture device to realize the same set of operations; 6.
Multiple types of machine tools and fixture devices are available. A single machine tool has the possibility to equip with different fixture devices, and vice versa; 7.
Three parameters for each type of machine tool are known: the mean time to failures (MTTF), the mean time to repair (MTTR), and the machine tool cost; 8.
Each fixture device enables access to multiple machining features and a datum system used to locate the part; 9.
A datum system has its machining features (used to locate the part) that have to be machined before the datum system is used; 10. The accessible machining features of each station decide the possible machining operations to be allocated to that station.

Production Information
The working time of the manufacturing unit (T A [hours/month]) can be obtained as: where T shift is the working hours per shift, N shift is the number of shifts per day, and D w is the working days per month. The expected demand interval is defined as [DM lb , DM ub ], where DM lb (DM ub ) is the lower (upper) bound of the expected demand per hour.

Manufacturing Information
Three kinds of feature are considered as follows. General features have two common machining constraints. One is "First roughing and then finishing", which means the precedence among the operations inside each feature is given. The other one is "First face milling and then drilling", which means that the drilling of the hole has to be completed after the machining of the one which the hole opens. For example, general features in the proposed case study can be classified into groups of features by the surfaces FGS = {FG 1 , . . . , FG i , . . . , FG SN }, where i = 1 . . . SN, and SN is the number of surfaces. For each FG i , i = 1 . . . SN, there are FNG i machining features with at least one operation that needs to be operated following its sequence. Operations in each FG i can be described as: In the equation, A = {a 1 , a 2 , . . . , a m } stands for the types of machining operations. In our example, {a 1 , . . . , a 3 } are possible machining operations of a surface (rough milling, semi-finish milling, and finish milling); while {a 4 , . . . , a 12 } are possible machining operations of a hole (drilling, gun drilling, core drilling, rough boring, semi-finish boring, finish boring, rough reaming, finish reaming, and tapping). Then, the element b k,j shows: b k,j = 1 a k have to be processed in feature j 0 a k does not exist in feature j The number of operations in each OFG i is: Features for each datum system (FDS i , with i = 1 . . . DSN) are features that have to be completed at the same station where the datum system has never been used before. According to ISO GPS and ASME GD&T standards, a datum system is a reference system for dimensional and geometric tolerances, and in machining it consists of features that are usually considered in relation to clamping and positioning the parts. For example, the datum systems of the cylinder block in this paper usually have one surface and two corresponding pin holes. DSN is the total number of datum systems. The operations of a datum system are combined as a single datum operation. So DSN is also the number of operations.
Special feature group (SFG i , i = 1 . . . SGN) is the set of features requiring special constraints, e.g., when operating two intersecting holes, both of them have to be machined at the same station and the longer hole should be machined before the shorter one; therefore, they are grouped into a special feature group, considered as a single special operation. SGN is the number of operations in the special feature groups.
Eventually, the number of machining operations of a part can be calculated as: Therefore, the station configuration group SCG = {SC 1 , . . . , SC i , . . . , SC SCN } is obtained. Each SC i shows a possible combination of F j and MT k . SCN stands for the total number of possible station configurations.

Relationship Information
Each feature group in the set FGG = {FG 1 , . . . , FG SN , SFG 1 , . . . , SFG SGN } has at least one station configuration to allocate the machining operations. Thus, the feature-station matrix FSM(i, j) is created to describe the relationship between FGG(i) and SCG(j), with i = 1..SN + SGN and j = 1 . . . SCN. It decides if operations from FGG(i) are accessible in the SCG(j): Whenever a station configuration is used, the features of its datum system must be machined at one of the previous stations in the line. Therefore, the datum-station matrix DSM(i, j) with i, j = 1 . . . SCN, is defined if operations from the datum features of SC i are accessible in the SC j : if datum features of SC i can be machined at SC j if datum features of SC i cannot be machined at SC j , i = j An extra situation happens at the first station, where the station configuration may consider the rough features (or pre-machined features) as the datum system. Therefore, DSM(i, i) is used to describe this situation as: if datum features of SC i do not need to be machined if datum features of SC i need to be machined (8)

Decision Variables
The considered decision variables are the following:

Relations and Constraints
To fulfil the real situation, the following relations and constraints should be verified.

Feature Group-Station Configuration Matrix
Feature Group-Station Configuration Matrix M fg−sc is mostly related to SI(s) with s = 1 . . . S. This matrix is used to show whether the machining operations of each feature group can be allocated to the selected stations:

Datum Feature-Station Configuration Matrix
Datum Feature-Station Configuration Matrix M df−sc with i = 1 . . . SCN is also related to SI(s) , with s = 1 . . . S. This matrix is used to show where the operations from each datum feature group can be allocated to the selected stations:

Station Constraints
Three station constraints are considered here. Firstly, the selected station configurations must have the ability to machine the operations from all the feature groups: Secondly, for the selected station configuration SI(k), with k = 2..S, there must be at least one pervious station that can machine the operations from its datum feature group: Thirdly, the first selected station configuration SI(1) must choose a fixture that can be directly used without machining the datum features:

Operation Constraints
Two main constraints are introduced for the operations. The first one is that each operation must be allocated to only one station: The second one is that each feature has its own sequence of machining operations. Thus, for any pair of machining operations of a feature o ∧ o = 1..OPN, respectively, realized in stations s ∧ s = 1..S, with o preceding o, the following constraints must be verified:

Buffer Capacity Constraints
The number of buffer capacities between stations should not exceed its limit, BA max , which ensures the limited spaces of each buffer area: 3.4.6. Investment Constraints C max stands for the maximum investment cost by the company. So, the total cost (C total ) must satisfy the maximum limit:

Objective Function
The two objective functions PR and C total are considered for the proposed problem. PR stands for the production rate of the machining line, the goal being its maximization. C total stands for the total cost that contains all the costs of machines and buffers, and the goal is its minimization. It can be calculated using the following equation, where MC is the cost of a machine tool, and BC is the cost of a single buffer element.
Tecnomatix ® Plant Simulation software developed by Siemens PLM Software is used to estimate the PR of solutions. As a discrete-event simulation software, its use has increased both in industry and recent research (see [30] as a review). The simulation software is integrated through the COM interface with the proposed algorithm developed in Matlab platform.
Therefore, the two objective functions considered are:

Solution Approach
In order to clearly explain the approach, we introduce the related logic diagram ( Figure 2). As shown in the figure, there are nine steps on the whole algorithm. Starting from "Definition of Basic Production Line", the initial information is prepared and defined. "Initial Population Creation" is the step used to randomly create the initial population that follows the given population parameters. The algorithm has a main loop whose ending policy is to check whether it has reached the maximum number of iteration times of the generations. During the main loop of the algorithm, individuals of the population in the current generation are decoded into feasible solutions through heuristic principles, and the two objectives (namely total investment cost and production rate) are then calculated for each solution. The total investment cost consists of machine tool and buffer capacity costs, while the production rate of each solution is estimated in a simulation environment. Hence, all the evaluated individuals in the current generation are sorted and ranked according to Pareto optimizing in NSGA-II [31]. Afterwards, some individuals from the worst solution set are eliminated to obtain the same population size and to form the final population of the current generation. Whenever the ending policy is not satisfied, the two NSGA-II operators (Crossover and Mutation) are introduced to generate the new population, which combines the current population with newly generated individuals. Therefore, the new population returns to the decoding step again, forming the main loop. After the ending policy is satisfied, the optimal results of the best solution set are finally presented.
The heuristic approach in this paper is more focused on the two core steps highlighted in grey in Figure 2, which are "Feasible Solutions by Decoding" and "Multi-objectives Calculation". Therefore, this section is organized as follows. Firstly, Section 4.1 presents the new logic diagram to show the complete logic of the proposed approach. Then, Sections 4.2-4.4 describe in detail the three main steps of the logic diagram (new heuristic principles, pre-screening strategy, and simulation strategy), respectively. Finally, in order to ensure a better performance of the algorithm, all the parameters are fine tuned in Section 4.5.
4.1 presents the new logic diagram to show the complete logic of the proposed approach. Then, Sections 4.2-4.4 describe in detail the three main steps of the logic diagram (new heuristic principles, pre-screening strategy, and simulation strategy), respectively. Finally, in order to ensure a better performance of the algorithm, all the parameters are fine tuned in Section 4.5.

New Approach Logic Diagram
The logic diagram of the heuristic parts (the dotted frame in Figure 2) can finally be summarized in Figure 3.

New Approach Logic Diagram
The logic diagram of the heuristic parts (the dotted frame in Figure 2) can finally be summarized in Figure 3.
Step 1: Definition of the current population set: Define the current population set Pop cur (M), where M = 1 . . . N ind , and N ind stands for the number of individuals of the population.
Step 2: Decoding of SI, OI and BA (Principles I-III): Decode the station configuration array SI, operation information array OI, and buffer allocation array BA, respectively, according to Principles I, II, and III from the original approach.
Step 3: Decoding CI (New Heuristic Principles): Decode the optimal configuration information array CI for the current individual through the new heuristic principles based on the demand interval in Section 4.2.
Step 4: Condition: Has the Population Set Pop cur (M) been fully explored? If so, go to Step 5; otherwise, M = M + 1, and go back to Step 2.
Step 5: Eliminating Solutions: Eliminate both the identical and worse solutions through the pre-screening strategy in Section 4.3.
Step 6: Multi-objectives Calculation with Simulation Strategy: Obtain the investment cost through Equation (6) and the throughput from the simulations based on the strategy as described in Section 4.4.
Step 7: Evaluated Solution Set: Obtain the evaluated solution set of current generation.  Step 1: Definition of the current population set: Define the current population set Pop (M) , where M=1…N , and N stands for the number of individuals of the population.
Step 2: Decoding of SI, OI and BA (Principles I-III): Decode the station configuration array SI, operation information array OI, and buffer allocation array BA, respectively, according to Principles I, II, and III from the original approach.
Step 3: Decoding CI (New Heuristic Principles): Decode the optimal configuration information array CI for the current individual through the new heuristic principles based on the demand interval in Section 4.2.
Step 4: Condition: Has the Population Set Pop (M) been fully explored? If so, go to Step 5; otherwise, M = M + 1, and go back to Step 2.
Step 5: Eliminating Solutions: Eliminate both the identical and worse solutions through the pre-screening strategy in Section 4.3.
Step 6: Multi-objectives Calculation with Simulation Strategy: Obtain the investment cost through Equation (6) and the throughput from the simulations based on the strategy as described in Section 4.4.
Step 7: Evaluated Solution Set: Obtain the evaluated solution set of current generation.

New Heuristic Principles
New decoding principles in this paper were built in order to meet the industrial requirements of the expected demand interval. The decision variables that make up the complete solution were grouped in the following four sets: configuration information (CI), station information (SI), operation information (OI), and buffer allocation (BA). The

New Heuristic Principles
New decoding principles in this paper were built in order to meet the industrial requirements of the expected demand interval. The decision variables that make up the complete solution were grouped in the following four sets: configuration information (CI), station information (SI), operation information (OI), and buffer allocation (BA). The structural encoding method was introduced where each individual can be illustrated as code P in Figure 4.  As shown in the figure, each code P consists of five vectors (P , P , P , P , and P ) and can be divided into three parts of information: station configuration, operation allocation and sequencing, and buffer allocation. During the first step of the decoding stage, the three parts are decoded into station configuration array ( SI ), operation information matrix (OI), and buffer allocation array (BA) through Principle I, Principle II, and Principle III based on the priority coding and constraints, where the three principles have ensured the feasibility of SI, OI, and BA. Firstly, SI is decoded by Principle I, where the three constraints in Equations (12)-(14) are satisfied. Then, OI is obtained by Principle II together with SI, where Equations (15) and (16)   As shown in the figure, each code P consists of five vectors (P 1 , P 2 , P 3 , P 4 , and P 5 ) and can be divided into three parts of information: station configuration, operation allocation and sequencing, and buffer allocation. During the first step of the decoding stage, the three parts are decoded into station configuration array (SI), operation information matrix (OI), and buffer allocation array (BA) through Principle I, Principle II, and Principle III based on the priority coding and constraints, where the three principles have ensured the feasibility of SI, OI, and BA. Firstly, SI is decoded by Principle I, where the three constraints in Equations (12)- (14) are satisfied. Then, OI is obtained by Principle II together with SI, where Equations (15) and (16) are fulfilled. Afterwards, BA is calculated by Principle III, satisfying the buffer constraint.
To handle the industrial requirement of the demand interval, new heuristic principles were proposed to decode the last element based on the demand interval and configuration information array (CI), which is the information used to decide the number of machine tools on each station of the production line. Considering the fact that the station configuration array (SI) has already been obtained, the number of stations (S), the fixturing device and part orientation, and the kind of machine tools of all stations are known. The investment cost of a single machine tool at station s is MC(SI(s)), where s = 1 . . . S, while the availability of the type of machine used at station s, A(SI(s)) can be obtained as the following equation Then, the feasible range of CI(s) for each station can be calculated in Equation (21), where CI(s) must be a positive integer: In a real situation, there is also the possible constraint of a minimum (CI min ) and maximum (CI max ) number of machine tools per station. Therefore, the final range of CI(s) for each station can be converted into Equation (22): , CI max , s = 1 . . . S (22) Since Equation (5) may provide multiple plans on the number of allocated machine tools at each station, the production line has more possible combinations for the configuration information CI. In order to further investigate these possible combinations for the current individual, all the possible CIs are generated. First of all, the total investment cost of each combination is calculated through Equation (23). Then, combinations that exceed the investment budget (Equation (24)) are all removed to fulfil the constraint of the maximum investment cost on the production line.
C total ≤ C max (24) Afterwards, in order to obtain the best of the remaining combinations, it is essential to find a proper function to evaluate them. As the considered problem is the combination of TLBP and BAP, we tried to introduce some common evaluating functions from TLBP to achieve the estimation of the solutions. Therefore, the two most typical objectives of line balancing were introduced: bottleneck time and deviation of allocation time.

1.
The bottleneck time of a production line is the maximum cycle time among all the stations. It represents the cycle time of the whole production line. However, in this situation, the solutions may have different levels of investment cost, which may make them incomparable simply through the bottleneck time; 2.
On the contrary, the deviation of allocation time has a more suitable mechanism in this situation. The deviation of allocation time stands for the deviation between the station cycle time and the mean cycle time at all stations of the production line. It evaluates the intrinsic quality of balancing for a single solution under the current investment cost of the manufacturing resources.
Therefore, the sum of the squares of the time deviation (SQ) was built here to evaluate all the combinations. As CTS(s), s = 1 . . . S of each combination can be calculated through Equation (20), the mean CTS(s) can be easily obtained as: In addition, due to the requirement of the demand interval, the solutions whose CT mean belong to the interval of [PT min ; PT max ] apparently have more advantages than the others. Therefore, the evaluating function of SQ of CTS(s) is proposed as: In Equation (9), there is a penalty mechanism that would enlarge the sum of squares of the deviations and provide such combinations with a worse evaluation. It happens for the combinations whose CT mean fulfils the two following conditions: either CT mean > PT max or CT mean < PT min .
Thus, the configuration information with the minimum value of SQ among all the combinations was considered as the final CI for the current individual. Finally, the individual was fully decoded.

Pre-Screening Strategy
As one kind of heuristic algorithm, NSGA-II generates the individuals more casually and evaluates all solutions of each generation. Although the heuristic principles jointly guarantee that each individual is translated into a feasible solution, the algorithm still has the chance to generate obviously worse solutions. As the majority of computing time of the whole algorithm is occupied with evaluating solutions through simulation, we decided to propose a pre-screening strategy to eliminate the mentioned solutions of each generation during the calculation. Thus, a pre-screening strategy here should have two abilities: 1 to identify the same solutions between two generations; 2 to roughly estimate the quality of each solution.
Firstly, since the coding system is priority-based, it may be not suitable to check if two individuals are identical or not directly based on their codes. So, the solutions from two neighbor generations will be compared after the decoding stage. The solutions are compared in pairs with their four key elements: configuration information array (CI), station information array (SI), operation information matrix (OI), and buffer allocation array (BA). To further reduce the computing effort, the comparison between OI matrix can be replaced by the array of total machining time at all stations MTM(s), where s = 1 . . . S.
The other proposal of the strategy was to roughly estimate the production rate or some other quality-related factors of the solutions before their evaluation through simulation, and then decide whether the solution is worth further evaluation through the simulation. Therefore, it was necessary to find a proper pre-screening objective function to evaluate the solutions. As discussed in the previous section, the deviation of allocation time provides the intrinsic quality of balancing that is independent of the manufacturing resources of the solution. Apparently, it is also suitable for pre-screening solutions with different levels of investment cost in each population. In addition, as SQ has already been obtained during the decoding stage, reusing it here would of course save the computational time of the algorithm. Therefore, the sum of the square of the time deviation (SQ) in Equation (9) was also used as the evaluation criteria to roughly estimate the quality of all individuals from each population.
After the evaluation of all solutions, a threshold function was introduced to make the decision on whether each solution should be eliminated before its simulation. In fact, the mechanism here is quite similar to the common optimizing direction when solving TLBP separately. The only difference is that TLBP uses this method to find the optimal solution and make the final decision, while we introduced it as a threshold to eliminate obviously worse solutions, thus avoiding the further search for BAP solutions in the invalid solution space.
Starting from the population in the initial generation (G 0 ), the default threshold (TR 0 ) is set as a maximum value that allows all or most solutions from the initial population (Pop 0 ) to have the chances to be evaluated through simulation. Then, the maximum value of SQ among the population Pop n of the current generation G n is used to define the threshold (TR n+1 ) for the population Pop n+1 of the next generation (G n+1 ). Therefore, the threshold can be calculated in Equation (27), where n stands for the generation, N gen is the total number of generations, and N plan is the population size:

Simulation Strategy
To enhance the stability on the throughput estimation of each solution, we performed the following experiment to study the uncertainty of the simulation and define the simulation strategy to use. To avoid a relevant computational effort of the proposed approach due to simulation, we considered three levels of the number of simulations (Sim_n) equal to 1, 3, and 5 times for each solution, and tested 60 randomly generated solutions. During the evaluation of each solution, the median value among multiple simulation results was used to evaluate the solution, which was the reason why we choose an odd number of times. For each solution, 10.000 replicas were run at the three levels.
We studied the three distributions of the 10.000 replicas for each of the 60 solutions, and present the detailed data in Appendix A. Apparently, as expected, the variability gradually reduced when the number of simulations on each solution increased, with the maximum range percentages equal to 2.2%, 1.5%, and 1.2% of the production rate, respectively, for Sim n = 1, Sim n = 3, and Sim_n = 5. The reduction in the range percentage going from Sim n = 1 to Sim n = 3 was much higher than going from Sim n = 3 to Sim n = 5. Noting that the additional simulation time required for the above two situations (Sim n = 1 to Sim n = 3 and Sim n = 3 to Sim n = 5) was the same, we decided to set Sim_n equal to 3 as the number of simulations among which to take the median value to estimate the production rate of each solution. Considering that the uncertainty could be considered symmetrically distributed, it was estimated to be ±0.75%.

Tuning of Algorithm Parameters
The parameters of NSGA-II are the size of the population N plan , the number of generations (N gen ), the percentage of crossover (P crossover ), the percentage of mutation (P mutation ), and the mutation rate (R mu ). To tune these parameters, we used a Taguchi approach considering the set coverage function C, presented by Zitzler et al. [32], as the response function. In order to compare two solution sets from Pareto optimization the set coverage was widely applied. For example, X and X were the two target sets, and the function C maps the ordered pair X , X into the interval [0, 1]: The value of C X , X equal to 1 shows that all solutions of X were dominated by or equal to solutions in X . The value of C X , X equal to 0 means the opposite situation where none of the solutions in X are dominated by or equal to solutions in X .
Using the set coverage, an L27 Taguchi design [33] was introduced to test the five factors with three levels (Low, Medium, High). Table 1 provides all the levels of tuning parameters, and Tables A1 and A2 in Appendix B show the detailed results of the experiment and the response table for means. Figure 5 provides the main effect plot for the means of the tuning results. Based on these results, the considered values of the parameters are shown in Table 2.

Industrial Case Study
This section presents the full case study to illustrate the design process of th manufacturing system for the four parts (engine blocks). In order to provide th possibility of future comparisons, the full data of the four parts used for optimization

Industrial Case Study
This section presents the full case study to illustrate the design process of the manufacturing system for the four parts (engine blocks). In order to provide the possibility of future comparisons, the full data of the four parts used for optimization can be found in this section and Appendices C-F. In this manufacturing system, there are four manufacturing units producing four kinds of engine blocks (Part A, Part B, Part C, Part D), which are shown in Figure 1. There are 49, 50, 57, and 58 different machining features with 84, 68, 109, and 96 different machining operations, respectively, for Part A, Part B, Part C, and Part D. In order to design and optimize the four manufacturing units of the system, their starting information was prepared according to the inputs listed in Section 3.2, which included "production information", "manufacturing information", "station information" and "relationship information".

• Production information
This plant is expected to work on a two-shift per day basis (8 h per shift), 300 days per year, so the annual working hours are 4.800. The cost of a single buffer capacity is CNY 10.000. Table 3 provides the production requirements of the four units in the manufacturing system for their target parts, including their total machining time and the expected demand range determined by the demand forecast of the company (DM lb is the demand lower bound, DM ub is the demand upper bound).

• Station information
In station information, machine tools, fixtures, and their datum system are included. Firstly, two kinds of 4-axis horizontal machining centers were considered for the four units and their parameters are shown in Table A3. There were 8, 9, 9, and 11 possible station configurations, respectively, for Part A, Part B, Part C, and Part D. Their synthetic representations are provided in Figure 6, and the detailed combination is shown in Table A4. be found in Appendix A, while the machining time of all the operations for the four parts are set out in Appendix B. •

Station information
In station information, machine tools, fixtures, and their datum system are included. Firstly, two kinds of 4-axis horizontal machining centers were considered for the four units and their parameters are shown in Table A3. There were 8, 9, 9, and 11 possible station configurations, respectively, for Part A, Part B, Part C, and Part D. Their synthetic representations are provided in Figure 6, and the detailed combination is shown in Table A4.  •

Relationship information
Feature-station matrix (row = feature group, column = station configuration), datum-station matrix (row = datum feature, column = station configuration) and special feature-station matrix (row = special feature, column = station configuration) are all built

• Relationship information
Feature-station matrix (row = feature group, column = station configuration), datumstation matrix (row = datum feature, column = station configuration) and special featurestation matrix (row = special feature, column = station configuration) are all built in Table A5, which illustrates the relationship and constraints between the feature groups and possible station configurations for the four parts.

Discussion
As we discussed in Section 2.3, referring to the industrial case study in [27], the simultaneous optimization of TLBP and BAP leads to a better overall solution (global optimization) with respect to the sequential optimization of the two problems (local optimization). To the best of our knowledge, reference [4] provides the only existing approach to the simultaneous optimization of the two problems. Therefore, we tried to compare this approach with the approach in [4] on the case study of four parts and the benefits of this approach are discussed in this section.
Based on the starting information in Section 5, the proposed algorithm with its tuned parameters was applied and compared with the approach in [4]. In order to ensure the robustness of the comparison, 30-time runs of both algorithms on each of the four engine blocks were considered. Therefore, 30 × 4 = 120 optimal sets (Pareto fronts) were obtained for each of the approaches.
Firstly, we compared these results at the level of optimal sets. Every optimal set from one algorithm was compared with the 30 optimal sets from the other one based on the set coverage function in Equation (28). Figure 7 reports the results. groups and possible station configurations for the four parts.

Discussion
As we discussed in Section 2.3, referring to the industrial case study in [27], the simultaneous optimization of TLBP and BAP leads to a better overall solution (global optimization) with respect to the sequential optimization of the two problems (local optimization). To the best of our knowledge, reference [4] provides the only existing approach to the simultaneous optimization of the two problems. Therefore, we tried to compare this approach with the approach in [4] on the case study of four parts and the benefits of this approach are discussed in this section.
Based on the starting information in Section 5, the proposed algorithm with its tuned parameters was applied and compared with the approach in [4]. In order to ensure the robustness of the comparison, 30-time runs of both algorithms on each of the four engine blocks were considered. Therefore, 30 × 4 = 120 optimal sets (Pareto fronts) were obtained for each of the approaches.
Firstly, we compared these results at the level of optimal sets. Every optimal set from one algorithm was compared with the 30 optimal sets from the other one based on the set coverage function in Equation (28). Figure 7 reports the results. As shown in the figure, the optimal sets of the new approach showed a significantly higher average value on their coverage over the ones from the original approach. Among the four cases, the new approach provided a quite high average value of the coverage on both Case A (78.71%) and Case C (83.27%), and a relatively high average value of the coverage on Case B (59.06%) and Case D (62.94%). On the contrary, the optimal sets of the original approach showed a very low average value on their coverage over the ones from the new approach (Case A 14.05%; Case B 27.04%; Case C 9.33%; Case D 22.01%). Therefore, the new algorithm presented in this paper showed great advantages on the quality level of the solution sets.
To further investigate the solutions inside the optimal set and the robustness of the algorithms, we performed an additional run on each of the four cases through both of the algorithms. The obtained optimal sets of the two algorithms applied on Case A are As shown in the figure, the optimal sets of the new approach showed a significantly higher average value on their coverage over the ones from the original approach. Among the four cases, the new approach provided a quite high average value of the coverage on both Case A (78.71%) and Case C (83.27%), and a relatively high average value of the coverage on Case B (59.06%) and Case D (62.94%). On the contrary, the optimal sets of the original approach showed a very low average value on their coverage over the ones from the new approach (Case A 14.05%; Case B 27.04%; Case C 9.33%; Case D 22.01%). Therefore, the new algorithm presented in this paper showed great advantages on the quality level of the solution sets.
To further investigate the solutions inside the optimal set and the robustness of the algorithms, we performed an additional run on each of the four cases through both of the algorithms. The obtained optimal sets of the two algorithms applied on Case A are shown in Figure 8 (Figures A2-A4 show similar results for Case B, Case C, and Case D, respectively). The horizontal and vertical axes of the figure represent the two objective functions: production rate and total investment cost, whose units were parts per hour and Million Yuan, respectively. Each single point represents a solution, where the red points are generated by the new algorithm, while the blue ones are obtained by the original algorithm: a point closer to the origin (left bottom) represents a better solution.
shown in Figure 8 (Figures A2-A4 show similar results for Case B, Case C, and Case D, respectively). The horizontal and vertical axes of the figure represent the two objective functions: production rate and total investment cost, whose units were parts per hour and Million Yuan, respectively. Each single point represents a solution, where the red points are generated by the new algorithm, while the blue ones are obtained by the original algorithm: a point closer to the origin (left bottom) represents a better solution. As shown in the figure: 1. Solutions in the optimal sets were distributed in groups and could be clustered.
Clusters mainly differed with respect to the number and type of machine tools, while the solutions inside a cluster mainly differed with respect to buffer capacity. Of course, the allocation of the operations also influenced the results; 2. In the comparison, the overall solutions from the new approach were better than those from the original approach in all four cases; 3. Considering the details of Case A, there were six clusters in the two optimal sets from both of the algorithms. The first two clusters with a lower production rate from the new algorithm showed slightly better solutions. Starting from the third cluster, the difference between the solutions from the two optimal sets was significantly larger. Considering the third clusters from both of the algorithms, the investment cost of the two clusters was quite similar and around 3.000 Million Yuan, while the production rates of solutions (around 9.9 parts/hour) from the new approach were relatively higher than those from the original approach (around 9.5 parts/hour). Starting from the fourth cluster, it is apparent that the solutions of the cluster from the original approach were fully dominated by those of the corresponding cluster from new approach. This behavior was similar in Case C ( Figure A3); 4. Considering the situation of Case B ( Figure A2) and Case D ( Figure A4), the improvement of the new algorithm with respect to the original one was less relevant.
To obtain a further understanding of these results, the optimal solutions of the additional run from both the approaches were clustered. The IP criterion with the investment cost per part was introduced and was calculated as: Therefore, the best solution (leader) for each cluster would be with the smallest IP value. The results of Case A are shown in Table 4 (new algorithm) and Table 5 (original As shown in the figure: 1.
Solutions in the optimal sets were distributed in groups and could be clustered. Clusters mainly differed with respect to the number and type of machine tools, while the solutions inside a cluster mainly differed with respect to buffer capacity. Of course, the allocation of the operations also influenced the results; 2.
In the comparison, the overall solutions from the new approach were better than those from the original approach in all four cases; 3.
Considering the details of Case A, there were six clusters in the two optimal sets from both of the algorithms. The first two clusters with a lower production rate from the new algorithm showed slightly better solutions. Starting from the third cluster, the difference between the solutions from the two optimal sets was significantly larger. Considering the third clusters from both of the algorithms, the investment cost of the two clusters was quite similar and around 3.000 Million Yuan, while the production rates of solutions (around 9.9 parts/hour) from the new approach were relatively higher than those from the original approach (around 9.5 parts/hour). Starting from the fourth cluster, it is apparent that the solutions of the cluster from the original approach were fully dominated by those of the corresponding cluster from new approach. This behavior was similar in Case C ( Figure A3); 4.
Considering the situation of Case B ( Figure A2) and Case D ( Figure A4), the improvement of the new algorithm with respect to the original one was less relevant.
To obtain a further understanding of these results, the optimal solutions of the additional run from both the approaches were clustered. The IP criterion with the investment cost per part was introduced and was calculated as: Therefore, the best solution (leader) for each cluster would be with the smallest IP value. The results of Case A are shown in Table 4 (new algorithm) and Table 5 (original algorithm), while the results of Case B, Case C, and Case D are shown in Tables A6-A11 in Appendix G. Moreover, an example of a detailed solution (the leader solution of cluster 5 in Table 4) is illustrated in Table A12. As shown in these tables, the results of Case A obtained by the new approach showed a much better spread of the production level (38.090, 43.008, 47.816, 52.719, 57.546, 62.334 parts per year) from the leading solution of the six clusters, and properly fulfilled the required interval of Case A (35.000, 60.000). Moreover, the best solutions from the new approach in Table 4 provided a quite similar and low investment cost per part IP (around 630 CNY/part), while those solutions from the original approach in Table 5 started from 633 CNY/part and gradually increased to 830 CNY/part. Apparently, the new approach proposed in this paper showed a better performance in Case A. Similarly, a better performance was obtained by this approach in Case C. The results from this approach provided a better spread of the production level (31.660, 35.209, 38.692, 42.273, 45.761, 49.028 parts per year) from the best solution of six clusters, meeting the required demand interval of Case C [28.000, 45.000]. The best solutions from this approach in Table A8 gave the lower investment cost per part IP (around 857 CNY/part), while the solutions from the original approach in Table A9 started from 858 CNY/part and gradually grew to 983 CNY/part. As for Case B and Case D, the two approaches showed a similar performance on the clusters with a lower production rate, while the new approach showed a better performance on the clusters with a higher production rate (especially the cluster with highest production rate). In addition, as shown in the tables, the new approach always provided a quite stable investment cost per part IP in all the four cases, which means that the new approach had a better searching ability within the solution space.
Therefore, the proposed approach showed great advantages and provided much better solutions in Case A and Case C, and also obtained better solutions in Case B and Case D.
In addition, to test the stability of both the algorithms, we ran 100 simulations for each solution of the 120 optimal sets and considered the median value for each solution as the reference point (namely PR med ). The percentage of deviation between the original production rate PR 0 and the reference point PR med was calculated as PR 0 −PR med PR med . The results of both the algorithms are shown in Figure A5. As shown in the figure, considering the 6934 solutions of the 120 optimal sets of the new approach, the median value of the deviation percentage was 0.58%, which was lower than the 0.89% obtained considering the 5678 solutions of the 120 optimal sets of the original approach. Thereforethis shows that the final solutions from the new approach were more stable than those from the original one.

Conclusions
This paper presented a heuristic approach for simultaneously solving transfer line balancing and buffer allocation problems considering the market demand uncertainty. We proposed to solve the overall problem using a combination of heuristic-based algorithm and NSGA-II. The new heuristic decoding principles enabled the industrial needs of the demand uncertainty to be met, which fulfilled the feasibility and improved the quality of each solution, thus achieving a better performance of the algorithm. During the decoding, the heuristic principles also provided each solution with a penalty mechanism to support the elimination of worse solutions. Based on this, we also proposed a pre-screening strategy to eliminate obviously worse solutions from each generation, thus greatly reducing the waste in computational time. Finally, we introduced and tested a simulation strategy to enhance the stability of the estimation of the throughput.
An industrial case study in a multi-unit manufacturing system producing four engine blocks was illustrated to both validate this approach and compare the novel approach with the original one. The comparison of the results showed that the novel approach provided significantly better results in Case A and Case C (with set coverage of 78.71% and 83.27%) and better performance in Case B and Case D (with set coverage of 59.06% and 62.94%), thus validating both the efficiency and effectiveness of this approach.
There are several aspects that might be further researched in the future. 1 From the problem point of view, future studies could consider an increased problem complexity to be closer to a real industrial situation, which would mean taking some other neighboring problems or more detailed information of the machining parameters or the machine tools into consideration. 2 From the solution point of view, the concept of multi-fidelity to evaluate BAP or other sub-problems might also be another relevant area for future research studies, which may enhance the efficiency of the methods.
Moreover, there are some suggestions that may help the potential industrial application in a similar situation. Simulations can provide much more information than the production rate used in this paper; thus, more relevant indicators such as the idle time of machine tools or the energy efficiency of the overall multi-unit manufacturing system could be introduced as objectives to evaluate the final performance of the manufacturing system according to the requirements of the decision makers.  number of simulations (Sim_n) equal to 1, 3, and 5 times for each solution. A total of 60 randomly generated solutions with 10.000 replicas per level were tested and are presented in the figure.
Chinese State Administration of Foreign Experts Affairs.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
This appendix provides the detailed simulation results considered in Section 4.4. Figure A1 includes summary reports for the percentage range of the three levels of the number of simulations (Sim_n) equal to 1, 3, and 5 times for each solution. A total of 60 randomly generated solutions with 10.000 replicas per level were tested and are presented in the figure.

Appendix B
This appendix provides the detailed data in relation to the tuning of the parameters of the proposed algorithm in Section 4.5.

Appendix B
This appendix provides the detailed data in relation to the tuning of the parameters of the proposed algorithm in Section 4.5. Table A1 describes the experiment results of L27 Taguchi of the proposed approach. Table A2 provides the response table for the means

Appendix C
This appendix provides the detailed information related to the station configuration of the four parts in the case study. Table A3 provides the parameters of the two machine tools considered in this paper, whose data were collected based on the previous production lines from the company. Table A4 provides the possible station configurations that could be used for all the four parts, whose combinations of machine type, datum system, part orientatio n, and accessible surface are shown in the table.

Appendix D
This appendix provides the detailed data in relation to the machining operations needed to produce the four parts (Part A, Part B, Part C, and Part D) considered in this paper. They include general operation groups (OFG), operation groups for datum system (FDS), and special feature groups (SFG).          Table A5 provides the detailed data in relation to the Feature Group-Station Configuration Matrix (FSM) for the four parts (Part A, Part B, Part C, and Part D) following Section 3.4.1. The station configuration can be found in Table A4 of Appendix C, while the information on the operation groups can be found in Appendix D.

Appendix F
This appendix provides the machining time of all the operations needed to produce the four engine blocks (Part A, Part B, Part C, Part D) considered in this paper.

Appendix G
This appendix provides the detailed comparison results discussed in Section 6. Figures A2-A4 show the comparison between the proposed approach and the previous approach at the level of optimal solution sets, respectively, corresponding to Cases B, C, and D, while Tables A6-A11 provide the comparison results at the level of solution clusters. Figure A5 shows the testing results on the stability of both the algorithms. Table A12 illustrates an example of the detailed leader solution of cluster 5 in Table 4

Appendix G
This appendix provides the detailed comparison results discussed in Section 6. Figures A2-A4 show the comparison between the proposed approach and the previous approach at the level of optimal solution sets, respectively, corresponding to Cases B, C, and D, while Tables A6-A11 provide the comparison results at the level of solution clusters. Figure A5 shows the testing results on the stability of both the algorithms. Table  A12 illustrates an example of the detailed leader solution of cluster 5 in Table 4 (Case A).

Appendix G
This appendix provides the detailed comparison results discussed in Section 6. Figures A2-A4 show the comparison between the proposed approach and the previous approach at the level of optimal solution sets, respectively, corresponding to Cases B, C, and D, while Tables A6-A11 provide the comparison results at the level of solution clusters. Figure A5 shows the testing results on the stability of both the algorithms. Table  A12 illustrates an example of the detailed leader solution of cluster 5 in Table 4 (Case A). Figure A2. Comparison for Case B: results from the optimal sets of the additional runs. Figure A3. Comparison for Case C: results from the optimal sets of the additional runs. Figure A3. Comparison for Case C: results from the optimal sets of the additional runs.       Figure A5. Deviation percentage distributions.