Multi-Objective Design Optimization of Flexible Manufacturing Systems Using Design of Simulation Experiments: A Comparative Study

: One of the basic components of Industry 4.0 is the design of a ﬂexible manufacturing system (FMS), which involves the choice of parameters to optimize its performance. Discrete event simulation (DES) models allow the user to understand the operation of dynamic and stochastic system performance and to support FMS diagnostics and design. In combination with DES models, optimization methods are often used to search for the optimal designs, which, above all, involve more than one objective function to be optimized simultaneously. These methods are called the multi-objective simulation–optimization (MOSO) method. Numerous MOSO methods have been developed in the literature, which spawned many proposed MOSO methods classiﬁcations. However, the performance of these methods is not guaranteed because there is an absence of comparative studies. Moreover, previous classiﬁcations have been focused on general MOSO methods and rarely related to the speciﬁc area of manufacturing design. For this reason, a new conceptual classiﬁcation of MOSO used in FMS design is proposed. After that, four MOSO methods are selected, according to this classiﬁcation, and compared through a detailed case study related to the FMS design problem. All of these methods studied are based on Design of Experiments (DoE). Two of them are metamodel-based approaches that integrate Goal Programming (GP) and Desirability Function (DF), respectively. The other two methods are not metamodel-based approaches, which integrate Gray Relational Analysis (GRA) and the VIKOR method, respectively. The comparative results show that the GP and VIKOR methods can result in better optimization than DF and GRA methods. Thus, the use of the simulation metamodel cannot prove its superiority in all situations.


Introduction
The fourth industrial revolution, known as industry 4.0, is considered the upcoming significant technology development as it allows customers to receive their products based on their expectations in terms of product varieties and quantiles [1]. Industry 4.0 can be attributed to its broadening focus on automation, decentralization, system integration, cyber-physical systems, etc. [2]. One of the basic components of Industry 4.0 is the Flexible Manufacturing System (FMS), which is an advanced production system that interconnects 1. According to the articulation of the preferences of the Decision Maker (DM). This first classification criterion is proposed by Rosen et al. [8]. Four groups of methods are possible and include the following: (1) a priori MOSO methods when the DM expresses their preferences before optimization is conducted; (2) a posteriori MOSO methods (in these methods, the DM selects a solution at the end of the search. Although this approach avoids the disadvantage of the a priori approach by taking into account preference information only at the end of the optimization process, it can lead to extremely high computational costs); (3) a progressive articulation of DM preferences (also named Interactive MOSO Methods) (the progressive approaches repeatedly solicit preference information from the DM to guide the optimization process). These methods enable DM to change his preferences during the optimization process by incorporating knowledge that only becomes available during the search. Interactive methods may be useful when simulation runs are expensive and the DM is readily available to provide input. Finally, the fourth group involves (4) non preference MOSO methods that operate without regard to the preference of DM.

2.
According to the research set and variables nature. This second classification criterion is proposed by Hunter et al. [7]. Three groups of methods are possible, including the following: (1) MOSO on finite sets, called Multi-Objective Ranking and Selection In the context of integer-ordered and continuous decision variables, we focus on methods that provably converge to a local efficient set under natural ordering. Furthermore, these methods of the three groups can also be viewed two groups according to the type of the final solution: global solution versus local solution [7]. The MORS methods provide a global solution, in which simulation replications are usually obtained from every point in the finite feasible set, and the estimated solution is the global estimated best. In addition, metaheuristics methods (named also random search) such as simulated annealing, Genetic Algorithms (GA), Tabu Search (TS), etc., also provide global solutions. Metaheuristics methods are efficient because they appropriately control stochastic error. However, the task is more challenging as it results in a number of solutions with different trade-offs among criteria, also known as Pareto optimal or efficient solutions. 3.
According to the use or non-use of metamodels. This third classification is proposed implicitly in many research studies such as in Barton and Meckesheimer [9], do Amaral et al. [10], etc. A metamodel or model of the simulation model simplifies the SO in two ways: The metamodel response is deterministic rather than stochastic, and the run times are generally much shorter than the original simulation. The metamodel is used to identify and estimate the relationship between the inputs and outputs of the simulation model, forming a mathematical function that is used to evaluate possible solutions in the optimization process. For example, Hassannayebi et al. [11] highlight that the adoption of metamodel-based SO in industry and service problems has grown due to its potential to reduce the number of simulation rounds necessary in the optimization process. Note that the MOSO methods, which are based on the metamodel, also provide a global solution such as that discussed in the second classification criterion.

FMS Design Literature Review
The study of Diaz et al. [12] presents a MOSO approach for a reconfigurable production lines subject to scalable capacities. The production line produces two product families and is composed of 18 workstations. The authors utilized a Non-Dominated Sorting Genetic Algorithm II (NSGA-II), a variant of GA to address the assignment of the tasks to workstations and buffer allocation for simultaneously maximizing the Throughput Rate (TR) and minimizing total buffer capacity.Červeňanská et al. [13] explored an MOSO of an FMS via a scalar simulation-based optimization method. The authors integrated a simulation with Design of Experiment (DoE) and Weighted Sum and Product multiobjective methods to optimize the total number of products, the Mean Flow Time (MFT), the Machine UTILization (MUTIL), and the average costs per unit of part. The modeled FMS produces two different products with eight workstations using parallel automated work machines.
The paper of Hussain and Ali [14] studied the impact of four design and control factors, control architectures, sequencing flexibility, buffer capacity, and scheduling rule on the performance of an FMS. The studied FMS is composed of six Computer Numerical Control (CNC) machines producing six different types of parts. The system is evaluated on the basis of make-span, average MUTIL, and the average Waiting Time (WT) of parts at the queue using the Taguchi-Grey multi-objective method. Apornak et al. [15] considered a multi-objective optimization of five performance measures in FMS. The authors addressed the optimal set of queues capacity, queues discipline, conveyor and transporter's speed, and operational setup times in an FMS with objectives of minimization of the average WT of raw materials, two average Process Times (PT), as well as the transporter and assembler product outputs. The studied FMS is composed of three work stations producing various kinds of seats for the freight cars. Using DoE, the authors simulated and collected the performance measure of 36 random scenarios. Regression analysis was then used to describe the metamodel of each performance measure. Consequently, the Response Surface Methodology (RSM) was applied to optimize the five objective functions.
Ahmadi et al. [16] proposed two Evolutionary Algorithms (EA): NSGA-II and NRGA are applied and compared to simultaneously combine the improvement of the make-span and stability of the schedule. This stability is evaluated by measuring the deviation of start and completion times of each job between prescheduled and realized schedule. The simulation is used to evaluate the state and condition of the machine breakdowns on a variety of manufacturing systems. Freitag and Hildebrandt [17] used a multi-objective simulation-based optimization to create a control strategy for an FMS by considering earliness and tardiness performance measures. This paper investigates the effect of 10 different attributes, which are the PT, the average PT of all waiting jobs, the Setup Time (ST), the average ST of all waiting jobs, the number of remaining operations, the time in system, the time in queue, the batch family size, the time until operational due date, and the average time until operational due date. The authors used the GA coupled with the simulation to solve the scheduling rule choice problem for a complex FMS.
Ammar et al. [18] investigated the size of the number of workers to be assigned to an FMS as well as the skills that each worker must have in a multi-objective optimization problem. The two objectives considered are minimizing the expected labor cost associated with the manufacturing team and minimizing the expected average task TR. The proposed multi-objective simulation optimization approach is applied to the design of teams of a manufacturing system; using the EA NSGA-II connected to a simulation model developed using Arena. Dengiz et al. [19] implemented a multi-objective optimization method of an FMS based on simulation through DoE, a regression meta-model, and the Goal Programming (GP) method. The authors have modeled and simulated by the ARENA simulation software an FMS with four workstations. Then, they applied the multi-objective optimization method to optimize the TR and MFT in the system by taking into consideration the number of operator, the velocity of material handling, the number of tool, scheduling rules, and the number of pallets as design and control parameters.
Using simulation results, Bouslah et al. [20] developed and solved a mathematical model based on RSM. The main objectives of the authors were to determine the optimal batch size, the optimal hedging level, and the economic sampling plan design, which minimized the average total holding cost, which includes the storage of the Work In Process (WIP) and final inventory stock, the average backlog cost, the average cost of sampling, the average costs of 100% inspection and rectification of the rejected batches, the average cost of transportation, and the average cost of replacement of non-conforming items sold to the consumer. However, the authors did not mention any details on the structure of the simulated manufacturing system. Iç et al. [21] considered a case study of simulation-based multi-objective optimization using the Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS) method hybridized with the Taguchi design technique. The studied production system is an FMS department composed of four CNC machining centers and producing three part types. The authors based their optimization case on the cycle time, TR, and work in queue as performance measures. In addition, they used five factors as decision variables. These factors are the number of cutting tools, the number of operators, the number of pallets, the velocity of transporter robots, and the pallet selection strategy.
The paper of Wang et al. [22] applies an MOSO method to a flexible shop scheduling problem. The two investigated objective functions are the minimum of the maximum PT and the minimum of the maximum machine load. The main considered constraints are the production resources and the technological process. The scheduling model of an FMS is established using simulation software and integrated to NSGA-II EA. In Zhang et al. [23], a hybrid method based on hybrid GA and TS is used to address a multi-objective FMS scheduling problem. Two objectives, which are the make-span and the starting time deviations, are considered to improve schedule efficiency and stability. A case of study of six machines FMS was studied with four different job arrivals rate and six different number of job arrivals.
Azadeh et al. [24] integrated simulations with the GP method and DoE technique to address a multi-objective scheduling problem of an FMS. The proposed method was applied on a real textile shop floor to minimize make-span and tardiness. The authors determined the decision parameters by using the DoE technique by estimating the effects the dyeing machine type, the temperature of the printing, the temperature and the number of center machines, and the scheduling rules through meta-modeling. Then, they used GP to find the optimal values of these decision variables, which are subject to a set of technical and managerial constraints. Um et al. [25] presented the simulation based multi-objective optimization of the design of an FMS with Automated Guided Vehicles (AGVs). Their principal objectives were to minimize congestion and utilization and to maximize TR based on many parameters including the number, velocity, and dispatch rule of AGV, part types, scheduling, and buffer sizes. In this paper, the authors considered a nonlinear programming method combined to evolution strategy. Nonlinear programming was used to determine the design parameters of the system through multi-factorial and regression analyses, and an evolution strategy was used to verify each parameter for simulation-based optimization.
Syberfeldt et al. [26] describe the use of Artificial Neural Networks (ANN) and EA as MOSO methods to the manufacturing cell at Volvo Aero. The two investigated objectives were the maximization of cell utilization and the minimization of overdue components considering the component inter-arrival times and due date as decision criteria. Kuo et al. [27] proposed a practical case of the Grey-based Taguchi method as a MOSO method for a company that provides integrated circuit packaging services. The authors aimed to optimize TR and cycle time performance for ink marking machines to avoid backlog of orders or lost customers, and the TR of the system must be increased. They based their methodology on five three-level control factors, which are the PT, the machine buffer size, the time between adjustment, the ratio of the adjusted PT to original PT, and the mean time between failures.
Oyarbide-Zubillaga et al. [28] focused on the determination of the optimal preventive maintenance frequencies for multi-equipment systems. The authors apply simulation and NSGA-II to the multi-objective optimization problem of preventive maintenance activities to minimize the system's cost and to maximize profit by considering the production speed, the percentage of unavailability of a machine due to corrective maintenance, and the fraction of time before and after the last maintenance as control factors. The system cost was defined as the sum of the preventive and corrective maintenance, the production speed lost, and the quality costs for each of the machines. Profit is the result of selling non-defective products. Park et al. [29] presented a method for determining the design and control parameters of an FMS with multi-objective performance via a fully factorial DoE, regression analysis and trade-off programming. A hypothetical FMS with six workstations was modeled and simulated. The number, speed, and dispatching rules of AGVs, in addition to the number of pallets, the buffer sizes, and the loading, routing scheduling rules, were considered as control parameters. These eight parameters were simultaneously determined by compromising performance measures of TR, delay, MUTIL, and WIP that are formulated using regression analysis.

The Proposed Conceptual Classification of MOSO for FMS Design
There are many MOSO methods applied for FMS design. According to the previous literature review, it is better to classify them in three main groups: Group A, Group B, and Group C, as detailed in Table 1. This classification is applicable regardless of the articulation of DM's preferences. It should be noted that all of the previous MOSO methods that are applied in the design of FMS are global solutions. A priori DM preferences are generally applied.

Group C No No
Iterative simulation and optimization using principally metaheuristics for random design research such as simulated annealing, genetic algorithms, etc. Only in this group, the articulation of the preferences of the DM is important. Table 2 summarizes the methods and techniques used in the MOSO methods used for FMS design. The presence of a cross "X" in a row and column intersection means that the research study stated in row use the method mentioned in column. It shows that all of the previous studies have applied a global solution method. These methods can be classified easily according the proposed classification in three groups (A, B, and C). Group C contains complex optimization techniques using metaheuristics, such as (GA, TS, EA, etc.). The performance of MOSO methods is not guaranteed because there is an absence of comparative studies. None of the previous studies has compared different MOSO methods.

The Objective of the Case Study
Our main contribution is to fill these gaps in the literature and to conduct a study of several relatively straightforward simulation-based FMS optimization methodologies that cover almost all categories of optimization methods classification. Our study investigates and compares the applicability and performances of the Goal Programming (GP), the Desirability Function (DF) method, the Grey Relational Analysis (GRA), and the VlseKriterijuska Optimizacija I Komoromisno Resenje (VIKOR) method.
All these methods are based on the DoE technique. They must be preceded by a design of experiments to program and sometimes analyze the simulation results. Moreover, these four multi-objective optimization methods have in common the type of preferences of DM; indeed, they are all based on an a priori decision of the DM for the choice of the objectives. On the other hand, the two methods GP and DF use the simulation-based metamodel technique and combine continuous and integer decision variables to solve the multiobjective optimization problem, while the two other methods, GRA and VIKOR, are based on the RS technique and use exclusively integer decision variables. The solutions reached by the GP and DF methods are then global and those reached by the GRA and VIKOR methods are local. In this study, we are interested in four multi-objective optimization methods in the context of FMS. An application on an FMS system will be used as a basis to compare the performances of these methods. It is mainly a matter of comparing the deviations between their results and the expected target values. Figure 1 describes in detail the adopted MOSO methodologies applied to FMS. These methodologies are essentially made up of three stages. Each of these stages consists of various steps. In the first stage, the primary step starts with FMS factors levels and performance measures selection and definition. Consequently, DoE is constructed, and the corresponding simulation models are developed using ARENA 14 discrete event simulation software. In the final step of this first stage, simulations are run to collect data for every studied performance measure. These simulation results are then analyzed in the second phase by one of the four adopted multi-objective optimization methods. The steps of this phase are discussed in detail in the following paragraphs. Finally, the optimum factors levels are adopted in the last stage of the multi-objective optimization method.

Materials and Methods
to the FMS of two consecutive parts, is generally generated by common probabilistic laws. In addition, parts are grouped into batches to reduce the machine's setup repetitions and transport times between work stations [31]. Furthermore, parts arriving at any work station are made to wait in a queue until the required machine becomes available. Once this required machine is idle, parts must be selected from the waiting queue based on scheduling rules [32][33][34]. As shown in Table 3, each of the considered FMS factors considered is studied with 2 levels.  The FMS investigated in this study is inspired by Pitchuka et al. [30]. An FMS is a manufacturing system characterized by a certain flexibility that allows the system to react in the case of changes. This flexibility is considered to fall into two categories. The first one, called routing flexibility, generally covers the system's ability to be changed to produce new product types. The second category is called machine flexibility, which consists of the ability to use various machines to perform the same manufacturing operation on a part.

•
To capture the FMS flexibility effect on its performance, this research adopts different machine LAYOUT (LAYOUT) for the studied FMS. Indeed, Functional Layout (FL) and Cellular Layout (CL) are the two most used machine layouts in FMSs. In FL, functionally similar machines are grouped into departments, and all machines of every department can perform production operations for any incoming part [31]. However, CL is made up of independent manufacturing cells. Each of these cells is made up of different machine types dedicated to the treatment of similar parts grouped into families. In addition, this work aimed also to measure the effect of part Batch Size (BS), part Inter-Arrival Time (IAT), and scheduling RULEs (RULE) on FMS performances ( Table 3). The IAT, defined as the difference between the arrival times to the FMS of two consecutive parts, is generally generated by common probabilistic laws. In addition, parts are grouped into batches to reduce the machine's setup repetitions and transport times between work stations [31]. Furthermore, parts arriving at any work station are made to wait in a queue until the required machine becomes available. Once this required machine is idle, parts must be selected from the waiting queue based on scheduling rules [32][33][34]. As shown in Table 3, each of the considered FMS factors considered is studied with 2 levels. • The FMS considered is composed of 8 machines grouped into 3 departments in FL and 2 cells in CL. The two departments "M" and "L" are composed of 3 machines each, while department "M" comprises only 2 machines. This MS is also characterized by two-part families composed of each of 2 part types. Each type of part requires 2 to 5 manufacturing operations (Table 4).
• The setup and processing times for each type of part are provided in Table 5. Setup times on every machine can be reduced or cancelled by the setup factor (δ) depending on the similarity of the successive parts family or type. Indeed, if successive parts belong to the same family, the subsequent part setup time must be reduced by a factor of δ = 0.5. On the other hand, if these successive parts have the same type, no machine setup is needed and the subsequent part setup time must be cancelled by a factor of δ = 0. Transfer times in the two layouts follow a statistical uniform law between 10 and 16 min. • To characterize the fluidity of parts flow in FMS, different optimization studies used WIP and MFT as major performance measures [31]. WIP has mainly been measured as the number of parts in the system, and MFT is simply obtained by averaging all durations between every part exit times and entry times in FMS. The TR of the production was adopted as the third performance measure. To evaluate TR, it is normal to measure the number of processed parts per unit of time. The maximization of such a measure of performance reflects the best use of material and human resources. To enhance the efficiency of FMS piloting, various optimization studies used the waiting and transfer times (WT and TT) as performance indicators, and they essentially aimed to minimize these two indicators.

The Simulation Model
FMS simulation models were built using Arena 14.0 software. The FL and CL models are composed of three parts: "Parts arriving", "Departments" or "Cells", and "System exit": • Parts enter to the system through a "Create" module named "Parts Arrival" in which the BS and IAT times are specified. Then, they are grouped into batches by a "Batch" module, named "Arrival Parts Grouping", to assign them their corresponding types through an "Assign" module named "Part Type". Due to the stochastic nature of their PT and ST, these batches are separated into unit products through a "Separate" module called "Parts Separation" to assign them each of their execution times through one of the four "Assign" modules named "Attribute Part i". However, a preliminary step must be performed through a "Decide" module called "Parts Sorting" to direct each type of product to the corresponding "Assign" module. The products then proceed through a "Batch" module named "Parts Grouping" before proceeding through the "Route" module named "Transfer to System" (Figure 2).
Machines 2022, 10, x FOR PEER REVIEW 10 of 26 durations between every part exit times and entry times in FMS. The TR of the production was adopted as the third performance measure. To evaluate TR, it is normal to measure the number of processed parts per unit of time. The maximization of such a measure of performance reflects the best use of material and human resources. To enhance the efficiency of FMS piloting, various optimization studies used the waiting and transfer times (WT and TT) as performance indicators, and they essentially aimed to minimize these two indicators.

The Simulation Model
FMS simulation models were built using Arena 14.0 software. The FL and CL models are composed of three parts: "Parts arriving", "Departments" or "Cells", and "System exit": • Parts enter to the system through a "Create" module named "Parts Arrival" in which the BS and IAT times are specified. Then, they are grouped into batches by a "Batch" module, named "Arrival Parts Grouping", to assign them their corresponding types through an "Assign" module named "Part Type". Due to the stochastic nature of their PT and ST, these batches are separated into unit products through a "Separate" module called "Parts Separation" to assign them each of their execution times through one of the four "Assign" modules named "Attribute Part i''. However, a preliminary step must be performed through a "Decide" module called "Parts Sorting" to direct each type of product to the corresponding "Assign" module. The products then proceed through a "Batch" module named "Parts Grouping" before proceeding through the "Route" module named "Transfer to System'' ( Figure 2). • As soon as one products batch arrives in one of the departments, it is separated into unit products and put on hold in the department queue via the "Hold" module named "Waiting Queue Department i". This queue is governed by a "Queue" module in which the scheduling rule must be specified. Once one of the department ma- • As soon as one products batch arrives in one of the departments, it is separated into unit products and put on hold in the department queue via the "Hold" module named "Waiting Queue Department i". This queue is governed by a "Queue" module in which the scheduling rule must be specified. Once one of the department machines becomes free, the selected waiting product is released from the "Hold" module. It then passes through a test, represented by the module "Decide" named "Machine Selection", which affects it toward this free machine. The processed products of the machines are grouped again in batches by the "Batch" module named "Grouping of Processed Parts Department i", which succeeds these machines. Finally, each batch of products is transferred to the next step in its production sequence through the module called "Route Department i" (Figure 3). • In the case of the CL simulation model, as soon as a batch of products arrives in on of the cells, it is directed to the first machine in its production sequence. This batch i then separated into unit products by a "Separate" module called "Separation Part Machine i". These products are then placed on hold in the queue of the machine via a "Hold" module named "Waiting Queue Machine i" until this machine become available. The choice of products from the machine queue is made according to th priority rule defined in the "Queue" module corresponding to this "Hold" module The processed products by one of the machines are grouped into batches via th "Batch" module called "Grouping Parts Machine i". This batch is transferred to th next machine in its production sequence via the "Intracellular Route Cell i" module By using this module, the transfer is performed in the cells, and the transfer time in this case is equal to zero. Each product with a completed production sequence mus be evacuated to the system's output section. Hence, the "Cell i Output Route" mod ule is used with a non-zero transfer time ( Figure 4). • In the two FL and CL simulation models, the machines are modeled by "Process" modules. In these modules, the transformation times are defined, which are function of PT and ST weighted by factor δ. Thus, Transformation time = PT + δxST. The valu In the case of the CL simulation model, as soon as a batch of products arrives in one of the cells, it is directed to the first machine in its production sequence. This batch is then separated into unit products by a "Separate" module called "Separation Parts Machine i". These products are then placed on hold in the queue of the machine via a "Hold" module named "Waiting Queue Machine i" until this machine becomes available. The choice of products from the machine queue is made according to the priority rule defined in the "Queue" module corresponding to this "Hold" module. The processed products by one of the machines are grouped into batches via the "Batch" module called "Grouping Parts Machine i". This batch is transferred to the next machine in its production sequence via the "Intracellular Route Cell i" module. By using this module, the transfer is performed in the cells, and the transfer time in this case is equal to zero. Each product with a completed production sequence must be evacuated to the system's output section. Hence, the "Cell i Output Route" module is used with a non-zero transfer time ( Figure 4). • In the case of the CL simulation model, as soon as a batch of products arrives i of the cells, it is directed to the first machine in its production sequence. This ba then separated into unit products by a "Separate" module called "Separation Machine i". These products are then placed on hold in the queue of the machi a "Hold" module named "Waiting Queue Machine i" until this machine bec available. The choice of products from the machine queue is made according priority rule defined in the "Queue" module corresponding to this "Hold" mo The processed products by one of the machines are grouped into batches v "Batch" module called "Grouping Parts Machine i". This batch is transferred next machine in its production sequence via the "Intracellular Route Cell i" mo By using this module, the transfer is performed in the cells, and the transfer ti this case is equal to zero. Each product with a completed production sequence be evacuated to the system's output section. Hence, the "Cell i Output Route" ule is used with a non-zero transfer time ( Figure 4).  • In the two FL and CL simulation models, the machines are modeled by "Process" modules. In these modules, the transformation times are defined, which are function of PT and ST weighted by factor δ. Thus, Transformation time = PT + δxST. The value of factor δ depends on the similarity of the types of products entering and leaving the machine. In fact, a module called "Selection Delta Value Machine Selection i" applies a test on all incoming products to the machine to look for the value of this factor. For this, it compares two variables named "Part Type" and "Part Family" defined in the two "Assign" modules, named "Part Type in Machine i" and "Part Type Out Machine i". If the two variables "Part Type" are identical, the module "Delta value machine selection i" directs the incoming product to the module "Assign" named "Delta Equal 0 Machine i" corresponding to the value of factor δ = 0. If the two variables "Part Type" are different but the two variables "Part Family" are identical, the module "Delta Value Machine i" directs the incoming product to the module "Assign" named "Delta Equal 0.5 Machine i" corresponding to the value of factor δ = 0.5. Otherwise, module "Selection Delta Value Selection Machine i" directs the incoming product to the module "Assign" named "Delta Equal 1 Machine i", corresponding to the value of factor δ = 1 ( Figure 5).
0, x FOR PEER REVIEW 12 of 26 incoming product to the module "Assign" named "Delta Equal 1 Machine i", corresponding to the value of factor δ = 1 ( Figure 5). • The leaving products batch proceeds through an "Assig" module, called "Output Performance Measures", for computing and updating all variables defined as performance measures. The acquired data are then stored in an Excel file using a "Readwrite" module for further treatment and analysis. Finally, the batches of products are evacuated from the simulation model via the "Dispose" module named "System Exit" (Figure 6).

Design of Experiments (DoE)
In this phase, we determine the number of distinct model settings to be run and the specific values of the factors for each of these simulation runs. There are many strategies for selecting the number of runs and the factor settings for each run include the following: random designs, combinatorial designs, sequential designs, factorial designs, etc.
Factorial designs are based on a grid, with each factor tested in combination with every level of every other factor. Factorial designs are attractive for three reasons: (i) The number of levels that are required for each factor is one greater than the highest-order power of that variable in the model, and the resulting design permits the estimation of coefficients for all cross-product terms; (ii) they are probably the most commonly used class of designs; and (iii) the resulting set of run conditions are easy to visualize graph-

•
The leaving products batch proceeds through an "Assig" module, called "Output Performance Measures", for computing and updating all variables defined as performance measures. The acquired data are then stored in an Excel file using a "Readwrite" module for further treatment and analysis. Finally, the batches of products are evacuated from the simulation model via the "Dispose" module named "System Exit" (Figure 6).
achines 2022, 10, x FOR PEER REVIEW incoming product to the module "Assign" named "Delta Equal 1 Machine sponding to the value of factor δ = 1 ( Figure 5). • The leaving products batch proceeds through an "Assig" module, called Performance Measures", for computing and updating all variables defined mance measures. The acquired data are then stored in an Excel file using write" module for further treatment and analysis. Finally, the batches of pr evacuated from the simulation model via the "Dispose" module named Exit" (Figure 6).

Design of Experiments (DoE)
In this phase, we determine the number of distinct model settings to be ru specific values of the factors for each of these simulation runs. There are many for selecting the number of runs and the factor settings for each run include the random designs, combinatorial designs, sequential designs, factorial designs, e Factorial designs are based on a grid, with each factor tested in combina every level of every other factor. Factorial designs are attractive for three reaso number of levels that are required for each factor is one greater than the hig power of that variable in the model, and the resulting design permits the est coefficients for all cross-product terms; (ii) they are probably the most comm

Design of Experiments (DoE)
In this phase, we determine the number of distinct model settings to be run and the specific values of the factors for each of these simulation runs. There are many strategies for selecting the number of runs and the factor settings for each run include the following: random designs, combinatorial designs, sequential designs, factorial designs, etc.
Factorial designs are based on a grid, with each factor tested in combination with every level of every other factor. Factorial designs are attractive for three reasons: (i) The number of levels that are required for each factor is one greater than the highest-order power of that variable in the model, and the resulting design permits the estimation of coefficients for all cross-product terms; (ii) they are probably the most commonly used class of designs; and (iii) the resulting set of run conditions are easy to visualize graphically for as many as nine factors [35].
The case study is about an FMS design with four factors, and each factor has two levels, as mentioned in the Table 3. Therefore, a 2 4 full factorial design was used to collect simulation results.

The GP Method
The GP is an optimization technique to solve problems with variety of objectives, which are generally incommensurable and often conflict each other in a decision-making horizon. The standard version of GP was first introduced by Charnes and Coper [31]. The GP model is based on an objective function formulated to find the most satisfactory solution that minimizes the total sum of positive and negative deviations from the level of attainment of the objectives levels (goals) set by the decision maker. This objective function is subject to physical and operating constraints of the system. The first type of constraints represents operating physical limits of the studied system. As for the second constraints, they are generally described by mathematical connections between the FMS factors and interactions and the performance measures to optimize. Hence, the principal purpose of the second stage of two first steps of the DoE-GP hybridization method is to build mathematical connections between FMS factors and responses. Statistical analyses are applied on the obtained simulation results to identify significant factors and interactions, and the relationships between the identified significant factors and interactions and the performance measures are translated to mathematical models by using the regression technique. In the third step of this stage, the GP model is developed setting the performance measures as goals and including other FMS constraints. Finally, this model is resolved using resolved using LINGO 18.0 software. The aim of this GP model is to find the most suitable levels of FMS factors that lower the total deviation of each performance measure from their respective target levels obtained in DoE.
The GP model takes the following form.
Moreover, it is subject to the following: ρx ≤ C (the operating physical constrain of the system) x j ≥ 0 (j = 1 . . . n), where the following is the case: 1. g i : The goal set for the ith objective for (i =1 . . . p) (the objectives here are the performance measures); 2.
x j : The jth decision variable for (j = 1 . . . n) (the decision variables here are the significant FMS factors and interactions); 3. a ij : The technological parameters (these parameters are the coefficients of the developed mathematical models relating the performance measures to significant FMS factors and interactions); 4.
ρ: The matrix of coefficients related to the physical FMS constraints; 5.
C: The vector of available physical FMS resources; 6.
δ + i , δ − i : The positive and negative deviations from the goals values.

The DF Method
The DF method is based on two steps. The first defines a desirability function by assigning values to responses that reflect their desirability. This involves transforming each value of the performance measure 'j' of experiment 'i', y ij , into a partial dimensionless desirability function d i , where 0 ≤ d i ≤ 1. This function includes the choices of the decision maker when constructing the optimization procedure.
A one-sided desirability transformation arises when the goal is to maximize or minimize the response, and two values A and B must be specified as the lower and upper limits. Equations (6) and (7) present the one-sided transformation equations that will be used for minimization and maximization goals, respectively.
The parameter ω j can be described as a power value or weight allocated according to the researcher subjective impression about the role of the response in the total desirability of the product.
A value of ω j equal to 1 implies that a linear desirability function is applied. If the value of ω j is less than 1, the obtained desirability function means that performance does not have to be close to the lower or upper limit, depending on the optimization goal, to have a higher desirability value. In contrast, if the value ω j is greater than 1, the desirability function implying that the performance has to be closest to the lower or upper limit, depending on the optimization goal, to have a higher desirability value.
To simultaneously optimize multiple performance sets, the individual desirability is combined using a geometric mean in the composite desirability.
A value of DF different from zero implies that all performances are in a desirable range simultaneously. In addition, a value of DF close to 1 means that the combination of the different criteria is globally optimal and the performances values are near the target values.

The GRA Method
Units of performance measurement are often different, so the influence of some of them may be neglected. This can also happen if some performance measures have a very wide range compared to others. In addition, if the expected optimization goals are contradictory, this will result in incorrect results in the analysis [36]. It is, therefore, necessary to normalize all performance values for each experiment in the first step of the multi-objective GRA-based optimization method's second stage.
In the developed DoE, for each of the "m" simulation experiments, "n" performance measures are measured. The ith experiment trial can be expressed as Y i = (y i1 , y i2 , . . . , y ij , . . . , y in ). Here, y ij is the value of the performance measure "j" of experiment "i". The term Y i can be translated into the comparability sequence X i = (x i1 ,x i2 , . . . ,x ij , . . . , x in ) using one of Equations (9) and (10), which are, respectively, used for larger-the-better and smaller-the-better objective values: x ij = y ij − y j y j − y j i = 1, 2, . . . , mj = 1, 2, . . . , n, x ij = y j − y ij y j − y j i = 1, 2, . . . , mj = 1, 2, . . . , n, where the following is the case.
After the normalization procedure, all x ij values, relative to the performance measures, will be scaled in [0, 1]. The Grey Relational Coefficient (GRC) is then computed to determine how close x ij is to x 0j = Max{x ij , i = 1, 2, . . . , m}. The larger the grey relational coefficient, the closer x ij and x 0j are. The grey relational coefficient can be calculated in the second step by the following: where the following is the case.
Once the entire GRC is computed, the Grey Relational Grade (GRG) is calculated in the third step based on the comparability and the reference sequence X i = (x i1 , x i2 , . . . , x ij , . . . , x in ) and X 0 = (x 01 , x 02 , . . . , x 0j , . . . , x 0n ) using the following: where ω j is the weight for the jth response, chosen by the decision makers. Of course, the sum of ω j is equal to 1.
In the final step of the GRA method, the GRG values are ranked in decreasing order. The optimal trial corresponds to the GRG maximum value.

The VIKOR Method
As in the case of the GRA method, which is based on GRG ranking, the VIKOR method is based on the computation of the VIKOR index and its ranking. In the first step of the VIKOR method, the ideal solution (A*) and the negative-ideal solution (A − ) are to be determinate. A* and A − represent, respectively, the maximum and minimum performance measure values of every experimental trial, and they are described as follows.
A * = Max y ij , i = 1, 2, . . . , m = y * 1 , y * 2 , . . . , y * j , . . . , y * n , A − = Min y ij , i = 1, 2, . . . , m = y − 1 , y − 2 , . . . , y − j , . . . , y − n , In the two following steps of the VIKOR method application the utility and the regret measures for the ith experimental trial, S i and R i respectively, are computed as follows: where ω j is the weight for the jth response, chosen by the decision makers. Of course, the sum of ω j is equal to 1.
In the fourth step, the VIKOR index of the ith experimental trial is computed as follows: where the following is the case.
Note that ν is the weight of the maximum group's utility. It is usually set to 0.5.
In the final step of the VIKOR method application, the VIKOR index values are ranked in decreasing order, and the optimal trials correspond to the maximum value.

Simulation Results
The case study is about an FMS design with four factors, each factor has two levels, as mentioned in Table 3. Therefore, a 2 4 full factorial design was used to collect simulation results. Each of the 16 simulation experiments was replicated 10 times. Simulation results show that a warm-up period of 10,000 min is needed, and models can then be run for 90,000 min. All final simulation results are provided in Appendix A. MFT simulation results are stated in Table A1, WIP simulation results are in Table A2, TR simulation results are in  Table A3, WT simulation results are in Table A4, and TT simulation results are in Table A5.

The GP Method
The use of GP as an MOSO method contains mainly four phases. The first phase is about the selection of the significant coefficient of the metamodel using Student's t-test. The second phase provides the final metamodel of each performance measure. The third and fourth phases concern the application of GP optimization.

1.
Determination of statistically significant FMS parameters: The main effects of the studied factors and interactions were analyzed in α = 0.05 of significance levels using the MINITAB statistical package (Table 6). Significant factors and interactions (p ≤ 0.05) are shown in bold.
With R 2 = 99.88% and R 2 (adj) = 99.87%, With R 2 = 97.86% and R 2 (adj) = 97.64%, Every constant in each of these equations corresponds to the average responses for each performance measure, and the coefficients assigned to the factors and interactions correspond to their respective effects.

3.
GP model formulation and resolution: We propose a GP model in which the selected performance measures are considered. The optimal configuration of decision variables minimizes the sum of penalties (dj). The parameter dj are deviations from the desired levels of the goals that are subject to series constraints. With the regression equations presented previously, the above-mentioned goal programming model can be stated as shown in Equations (32)-(40): which are subject to the following. −61.61 LAYOUT and RULE are binary (1 or 2), The objective G MFT , G WIP , G TR , G TT , and G WT goal values were fixed basing on the experimental design results.

4.
The GP model was solved using the mathematical software LINGO 18.0. The best value of the objective function was found to be equal to 136.99 and was obtained for the following levels of the studied factors: LAYOUT = CL, RULE = FCFS, IAT = 25, and BS = 5.

The DF Method
Applying Equations (6) and (7) for the studied performance measures, the individual desirability functions 'd' are very close to 1.0, as shown in Figure 7. Furthermore, Figure 7 illustrates the effect of each factor (columns) on the FMSs' performance measures and the desirability of the composite (rows). The red vertical lines and the corresponding numbers in red indicate the levels of optimal factors. The blue horizontal lines and the corresponding numbers in blue represent the values of the performance measures corresponding to the levels of optimal factors. Each of the performance measures is accompanied by the corresponding desirability function values 'd i '. In addition, the first row provides the value of the composite desirability 'DF', as presented in Equation (8), corresponding to the levels of the optimal factors. The obtained DF is equal to 0.984, which represents an ideal case of optimization. To obtain this desirability, the factors' levels must be set to the values shown below the global solution in Figure 7. That is, BS = 5, IAT = 5, LAYOUT = CL, and RULE = SPT.

The GRA Method
Based on Equations (9)-(17), the simulation results were normalized, and the GRC and GRG were calculated (Table 7). Once GRG was ranked, it appears that the optimum performance measures were obtained for the factor levels LAYOUT = CL, RULE = SPT, IAT = 5, and BS = 5. The row in bold in Table 7 indicate the optimal solution obtained using the GRA method, which has a rank equal to 1.

Discussion
The application of the four optimization methods in the context of FMS shows good results for four of the five performances in the case of the two methods GP and VIKOR, and only two performance measures for DF and GRA methods (Table 9).

The VIKOR Method
Based on Equations (18)- (26), the utility and the regret measures as well as the VIKOR index were computed (Table 8). Once the VIKOR index was ranked, it appears that the optimum performance measures were obtained for factor levels LAYOUT = CL, RULE = SPT, IAT = 25, and BS = 5. The row in bold in Table 8 indicates the optimal solution obtained using VIKOR method, which has a rank equal to 1.

Discussion
The application of the four optimization methods in the context of FMS shows good results for four of the five performances in the case of the two methods GP and VIKOR, and only two performance measures for DF and GRA methods ( Table 9).
The results show that the MFT, WIP, TT, and WT performance measures met their targets for GP and VIKOR methods. Indeed, they all show relatively minor deviations from their target values. Only the deviation of RT reaches, respectively, −63.53% and −81.87% in the case of these two methods. On the other hand, in the case of DF and GRA methods, only the optimal values of TT and TR were close to their corresponding targets, while the deviations between the achieved values and the objective values in the case of WIP, WT, and MFT can reach +242.529%, +107.377%, and +83.293% respectively. Hence, the optimization results can be considered satisfactory in the case of GP and VIKOR methods. Meanwhile, it was not the case for the optimization results of the DF and GRA methods. The two methods GP and DF require a higher level of analysis effort than the two methods of GRA and VIKOR. Indeed, in addition to the modeling and development of the simulation models, which is a common point to the four compared optimization methods, as well as the planning of experiments with the DoE method, the two methods GP and DF require relatively higher levels of expertise in the use of the analysis software MINITAB and LINGO. On the opposite side, the two methods GRA and VIKOR only need the development of the equations on Excel, which is within the reach of the majority of DMs. This has an impact on the applicability of the MOSO method. Table 10 summarizes the performance of the four MSOSO methods being compared in this study. Signs '+' and "−" are assigned to the optimization methods based on their achieved optimization results and their applicability. A '+' is assigned to each method resulting in a good optimization result, which is expressed by reasonable or small deviations. On the other hand, a "−" is assigned to each method that leads to an optimization result characterized by high deviations. For applicability, a "−" is assigned to each method that requires a high level of analysis and expertise. In the opposite case, a '+' is assigned to this optimization method. These methods are then classified according to assigned signs. Any method obtaining two "+" signs will be considered the most efficient. On the other hand, if it obtains two "−" signs, it will be considered as the most mediocre one. In the case where the optimization method obtains both signs "+" and "−", the classification gives priority to the obtained optimization result. The best method is VIKOR, which belongs to group B in the proposed classification. It is followed by the GP method, from group A, since it reaches good optimization results, although it requires a considerable analysis effort. The GRA method, from group B, comes in third rank and the DF method, from the group A, closes the classification at the last rank. This classification shows that the use of optimization methods based on a metamodel does not always produce the best results.

Conclusions
Various MOSO methods have been presented, developed, and used in the literature. These methods have been the subject of numerous classifications. However, the performance of these methods is not guaranteed due to the lack of comparative studies. Moreover, these classifications have been very diverse and are rarely related to the specific domain of manufacturing systems.
The objective of this research is two-fold. First, we proposed a new conceptual classification of MOSO methods applied to the context of MFS design. Second, four MOSO methods are selected according to this classification and compared through a case study related to the FMS design problem inspired by the literature. This comparison is based on the quality of the optimal solutions obtained by these methods as well as the degree of difficulty of their applicability through the necessary analysis effort and the degree of expertise of the user of these methods. All these studied methods are based on DoE. Two of them are metamodel-based approaches that incorporate the GP and the DF, respectively. The other two methods are not metamodel-based approaches and incorporate GRA and VIKOR, respectively. The comparative results show that the VIKOR method can result in a better optimization than GP, GRA, and DF methods in that order. It is clear, thus, that the use of MOSOs based on meta-models does not produce the best solution in all situations.
This research compares four MOSO methods applied in the context of FMS design. Some future research perspectives should be addressed:

•
In this study, four MOSO methods are compared. Two methods belong to group A of the proposed new classification, while the other two belong to group B. The extension of the current comparison to other MOSO methods belonging to group C is the first objective of our interesting perspectives.

•
The studied MOSO methods have been applied on a model of an FMS inspired from the literature. This model has six machines grouped in two cells in the CL and three departments in FL. In addition, this FMS processes only four products grouped into two families. Extending the comparison performed in this study to real and more complex FMSs to evaluate the reliability of MOSO methods is the second objective of our interesting perspectives. The application of the compared MOSO methods proceeds through different steps to generate optimization solutions. These steps usually require the intervention of a user to transfer the results from one step to another. The integration of these analysis and optimization steps into the simulation software, as in the case of the OptQuest tool in several simulation tools, would be a very interesting perspective.