Improved multiobjective differential evolution with spherical pruning algorithm for optimizing 3D printing technology parametrization process

Multiobjective optimization approaches have allowed the improvement of technical features in industrial processes, focusing on more accurate approaches for solving complex engineering problems and support decision-making. This paper proposes a hybrid approach to optimize the 3D printing technology parameters, integrating the design of experiments and multiobjective optimization methods, as an alternative to classical parametrization design used in machining processes. Alongside the approach, a multiobjective differential evolution with uniform spherical pruning (usp-MODE) algorithm is proposed to serve as an optimization tool. The parametrization design problem considered in this research has the following three objectives: to minimize both surface roughness and dimensional accuracy while maximizing the mechanical resistance of the prototype. A benchmark with non-dominated sorting genetic algorithm II (NSGA-II) and with the classical sp-MODE is used to evaluate the performance of the proposed algorithm. With the increasing complexity of engineering problems and advances in 3D printing technology, this study demonstrates the applicability of the proposed hybrid approach, finding optimal combinations for the machining process among conflicting objectives regardless of the number of decision variables and goals involved. To measure the performance and to compare the results of metaheuristics used in this study, three Pareto comparison metrics have been utilized to evaluate both the convergence and diversity of the obtained Pareto approximations for each algorithm: hyper-volume (H), g-Indicator (G), and inverted generational distance. To all of them, ups-MODE outperformed, with significant figures, the results reached by NSGA-II and sp-MODE algorithms.


Introduction
The product development conception during the manufacturing systems needs optimal experiment design techniques which can help to systematically develop experiments by incorporating important information to estimate parameters with higher accuracy. Normally, to optimize parameters related to additive manufacturing, the design of experiments methods has been used. However, when problems become complex and several objectives are taken, the optimization method allows obtaining faster and trustworthy solutions. The solution to a multiobjective optimization problem is normally not unique since the best solution for all objectives does not exist. There is a set of good solutions referred to as non-dominated solutions (none is better for all objectives) that define the Pareto set and the Pareto front (objective values for the Pareto set solutions) . Therefore, it is important to combine analytical and heuristic techniques for solving complex problems such as technical parameters in the manufacturing processes (Tervo et al., 2003). For decision-making, it is necessary to obtain the current reality, relevant metrics, and optimized measures for the best decision (Canciglieri et al., 2015;Chen & Zhao, 2016).
In this case, statistical techniques can be assumed for the design of experiments (DOE), as optimal experiment design techniques can help to systematically develop experiments by incorporating important information to estimate parameters with higher accuracy (Chen & Zhao, 2016). Based on market competitiveness, quality, and expectations of customers, DOE has been used to investigate hideaway causes of process variation. This strategy allows mapping the effects of hidden variables as well as to study the impacts of such effects during both process design and development. In a few words, DOE brings the range of experiments from uncontrollable factors that are introduced randomly to carefully control these factors (Antony et al., 2014). Wohlgemuth et al. (2020) integrated the statistical design and mathematical programming approaches to improve the productive efficiency of a set of decision-making units for logistic operators. Tervo et al. (2003) argued that, when a parameter appears nonlinearly in the model, it can be considered that an optimally designed experiment depends on the current estimated value. Wiecek et al. (2016) asserted that multiobjective programming allows optimizing decision problems through multiple objective functions. Ehrgott et al. (2018) argued that the task of optimizing is regularly used to define a robust model that gives the best solution for modeling and solving optimization problems with uncertain data. In this way, in a robust optimization process, both evolutionary and population stochastic optimization techniques can be applied to increase accuracy, convergence and at the same time as reducing the computational cost to obtain satisfactory results in terms of the decision variables (Consigli et al., 2020;Wang et al., 2020). Rao et al. (2017) assert that the nondominated sorting genetic algorithm II (NSGA-II) is one of the most popular multiobjective optimization algorithm used for optimizing machining operations parameters in minimizing or maximizing the machining performances. Bhavsar et al. (2015) solved a micro-milling of cemented carbide problem using multiobjective optimization NSGA-II algorithm to optimize the surface roughness and material removal rate produced. Alvarado-Iniesta et al. (2017) integrated the NSGA-II algorithm with machine learning and multicriteria decision-making techniques to optimize the dimension of the plastic product, cycle time, and packing pressure in the injection molding process. Alvarado-Iniesta et al. (2019), solved a many-objectives design problem associated with plastic gear using an NSGA-III algorithm to optimize seven objectives in a plastic injection molding process.
According to Salomon (1998), gradient-based algorithms can be applied only to continuously differentiable objective functions. But if either the objective function is not continuously differentiable or if the function is not (completely) given due to limited knowledge, which often occurs in real-world applications, the designer must resort to other methods, such as evolutionary or population-based algorithms.
Applications of population-based algorithms can be found in the literature in the field of DOE problems. Shyu et al. (2004) introduced a meta-heuristic based upon the Ant Colony Optimization (ACO) approach, to find approximate solutions to the minimum weight vertex cover problem, where the objective is to find a vertex subset whose total weight is minimum subject to the premise that the selected vertices cover all edges in the graph. Shih et al. (2014) presented a framework that adds a variable selection step prior to run a computer experiment-based optimization considering data mining methods. Additionally, using principal components analysis and multiple testing based on the false discovery rate for variable selection, the proposed method, which was dedicated to DOE, was applied to an airline fleet assignment case study. The Multiobjective Particle Swarm Optimization (MOPSO) algorithm was used to optimize the relation between sustainability practices and performance measurement of small and medium enterprises of plastic injection (Abdelaziz et al., 2018). El-Hajj et al. (2020) solved a vehicle routing problem using a MOPSO algorithm to maximize the total amount of profit, and the travel time limit is set. Trivedi et al. (2020) presented a review of recent PSO variants and solved benchmark and industrial optimization problems using a simplified PSO algorithm.
When it comes to evolutionary algorithms, Differential Evolution (DE) became one of the leading metaheuristics. DE consists of methods that operate over survival-of-the-fittest principles (Storn & Price, 1997). As an optimization tool, DE has presented its extension to single and multiobjective problems, emerging the term Multiobjective Differential Evolution (MODE). Zhang et al. (2013) presented an approach of MODE applied to DOE, where the combination of those techniques provided insights into the optimization of ferrite magnet machine. Validi et al. (2020) integrated DOE with a bi-objective metaheuristic to optimize in two phases a sustainable model that minimizes CO 2 emissions from transportation and total costs incurred in facilities and the transportation channels. Stokes et al. (2020) presented MODE applied to DOE as guidance for practitioners to find different types of optimal designs for various statistical models in chemistry applications. Reynoso-Meza et al. (2010) proposed the spherical pruning (sp) mechanism and adapted it to the classical MODE to ensure diversity of Pareto solutions, creating the sp-MODE. The new algorithm was subsequently applied to problems involving PID-controller tuning, considering the optimization of many-objectives (+ 4) simultaneously (Reynoso-Meza et al., 2010). The sp-MODE has also made its way into areas of engineering, being used by Hamdy et al. (2016) in a benchmark comparison of six different algorithms, and how well they perform over nearly-zero-energy-building design problems, as well as Camilotti and Freire (2020), utilized it alongside the NSGA-II in a multi-objective building design optimization problem. In both studies, the capacity of the sp-MODE in obtaining well-distributed and diverse Pareto solutions is accentuated.
This research proposes a hybrid approach to optimize the parameters, integrating the DOE, and multiobjective optimization methods. It has been developed a hybrid approach as an alternative to classical parametrization design used in the 3D printing process, by combining it with a proposed multiobjective differential evolution with uniform spherical pruning (usp-MODE) algorithm. The proposed modification in the algorithm has the objective of creating a more uniform sphere surface grid to be used in the spherical pruning mechanism. This work also presents an integrated approach to select the Pareto front inside a range of optimal solutions considering boundaries to the objective functions of the 3D printing process. The main process parameters and their influence on surface roughness, dimensional accuracy, and mechanical resistance of a 3D printing prototype specimen were analyzed.

Theoretical background
This section presents a brief literature review associating the machining process, design of experiments, multiobjective optimization, and fundamentals to reach the proposed usp-MODE algorithm. The main idea of this section is to provide a reasonable overview of the integration of these approaches.

Additive manufacturing and machining processes
In 1981, additive manufacturing (AM) was introduced regarding a solid printed model in Nagoya Municipal Industrial Research Institute. A few years later, stereolithography, which is the technique of creating models by curing a liquid photopolymer resin using UV lasers was presented. The previously mentioned technology, which was proposed by Charles W. Hull, was commercialized as the first rapid prototyping system (Hull, 1984). At the beginning of the 90's, Scott Crump invented the Fused Deposition Modeling (FDP), which is the basis of the current personal 3D printers (Crump, 1992). Medical applications using 3D printing started at the end of the 90's, while in the early 2000's the first commercially viable Selective Laser Sintering (SLS) machine, and the possibility of multiple materials mixtures were introduced (Beaman et al., 2020). Today, research-based new materials for 3D FDP printing considering injection technology were focused on the avoidance of hazardous resin materials. Moreover, recent developments have allowed laser-cured resin processes to emerge onto the commercial and home-user market.
The 3D printing technology was incorporated into the industrial environment and has been used in production systems due to the following strengths: the ability to produce complex and detailed three-dimensional shapes; reduced lead time for unique parts; installation possibility in non-industrial environments; low cost; and general savings compared to common manufacturing systems. The machining process has extended concepts of printing by giving the chance to build physical objects with any geometry in real or scale sizes just based on drawings or models designed in a CAD tool, as presented in (Canciglieri et al., 2015).
The concept of the additive manufacturing process in the industry came from the idea of printing final versions of finalized products. Based on the injection process to create layers that are deposited on a platform that moves in a plane, which thickness is given by a third axis, a laser pointer is guided along with the resin by moving mirrors. The resolution and accuracy to manufacture products today are so high that it becomes possible to construct highly complex parts that are no larger than a grain of sand. To reach the desirable efficiency, regarding the trade-off between time and quality of 3D printing products, the problem associated with the proper configuration of the machine became an issue discussed in academic and industrial environments.
The importance to build complex structures using the machining process and the potential of modeling by finite elements method to predict droplet size for manufacturing were already discussed and could be found in the specialized literature (Fernandes et al., 2017).
Due to the flexibility of using several materials to build 3D elements, the additive manufacturing research field faced new challenges, as machining configuration focusing on the optimization of the part being printed, and the optimization of the process in terms of time and scheduling inside the industry.

Design of experiment
Until the 70 s, experimental studies were constructed considering combinatorial rules. The classic DOE was developed by Sir Ronald A. Fisher in 1920, in 1970 experiments using optimization algorithms started to appear (Antony, 2014). The factorial design is indicated for studying the effects of two or more controllable factors (input variables) in the output variables of the process. The usage of meta-heuristic for optimization of experimental designs aims to improve quality or to optimize a model into the solution space using a convergence of the algorithm with accuracy increased and reduced costs (Validi et al., 2020).
In the design of experiments, a finite number of designs in the design space are simulated using prescribed settings of the design variables and system evaluation routines. It is used for process improvement, for minimizing the differences between results and nominal specifications, and for reducing both variability and overall manufacturing costs. In many manufacturing systems, this technique can be applied to select the best process parameters to increase robustness, evaluate alternative materials, and reduce the variation that affects the performance of the process or product (Sant'Anna, 2015). A robust process can be defined as a process that provides a reasonable response to input variables regardless of the external noises and with minimal resources to ensure a reliable product. The robustness of the process occurs due to a good experimental design that rigorous controls for materials composition (Zhang, 2009).
Regarding the number of factors, there are many different approaches for DOE analysis. The 2 k factorial design used in this research seeks to determine the optimum levels of the process for conducting several factors-at-a-time, explaining the significant effects in the process. This factorial design adopts k input variables with two levels each, one high and one low, and these levels may be quantitative or qualitative. These experiments are an important class of planning because the number of trials involved in conducting these experiments is relatively small (Dumas et al., 2014).
The main objective of applying this technique is to answer questions regarding the general behavior of the output variables within the range of input variables and to map regions of high performance. For this, three steps are important: (i) to plan the experiments by distributing accordingly the runs; (ii) to estimate the coefficients of the output variable; and finally, (iii) to explore the output variable region, finding the optimal level of the input variables that maximizes or minimizes the output.

Multiobjective optimization
In a generalized view, optimization can be interpreted as the search for the most appropriate solution for a given problem, given that the solution obeys a set of given constraints. From , basic concepts of optimization might also be viewed as the search for the set of values that minimizes a given function (Definition 1).

called a global minimum if and only if
In engineering, these problems can quickly reach high levels of complexity, considering the number of design variables considering and the underlying dynamic of the system being optimized (Rao, 2009). Moreover, multiple criteria are usually taken into consideration during a real-world problem, with them being intrinsically conflicting with each other (Definition 2).

Definition 2 (General Multiobjective Optimization Problem (MOP)): Find the vector.
T which will satisfy the m inequality constraints the p equality constraints and will optimize the vector function

Definition 3 (Pareto Optimality):
The conflicting nature between optimization objectives transforms the unique solution that the problem would otherwise have, in a set of non-dominated (Definition 4) solutions with different compromises for each criterion.

Definition 4 (Pareto Dominance):
The goal of a MOP is usually to find the optimal set (Definition 5) of non-dominated solutions, that cannot be improved any further by any other feasible solution. While the domain of the solutions forms the Pareto Set (Definition 5), their image on the objective space represents a curve denominated Pareto Front (Definition 6).

Definition 5 (Pareto Set) : For a given MOP
Definition 6 (Pareto Front) : For a given MOP ⃗ f (x) and Pareto Optimal set P * , the Pareto Front (PF * ) is defined as

Evolutionary algorithms and differential evolution (DE)
One of the preferred ways to tackle such problems is using evolution-based meta-heuristics, a practice that has gained popularity in the last years (Canellidis et al., 2016;Dumas et al., 2014;Eiben & Smith, 2015). Algorithms in this category provide a good balance between exploration and exploitation of the search space, abstracting themselves from problem-specific traits, and performing well over problems that classic methods would have trouble solving, if at all (Talbi, 2009). During the last decades, many metaheuristics have been developed, with the Differential Evolution (DE) (Storn & Price, 1997) being one of them.
In a similar manner to other evolutionary algorithms, the DE treats the input vector of design variables as an individual in a population of solutions, while their respective value over the objective function represents their fitness. Both mutation, crossover, and selection operators are an integral part of the algorithm, used as means to generate new individuals as well as preventing the gene pool from becoming stagnant (Storn & Price, 1997). At every iteration, every solution in the population acts as a parent x p , while three other random ones are used to create a mutant solution x m . Just as the algorithm's name implies, the differential term between the second and third solutions is multiplied by the differential weight F, and subsequently added to the first one (Eq. (9)).
The mutant and parent solutions are then combined through the crossover operator and a crossover chance Cr , resulting in a new child solution, as shown in Eq. (10). Finally, once the child solution x c is constructed, the selection operator is applied to both it and the parent solution, where the best fit one, under the objective function, is selected to continue to the next generation. The process is repeated until the desired stop criterion set by the user is reached, such as the number of generations or function evaluations (Storn & Price, 1997).
Due to its simplicity, many variants of the DE were proposed in these past years, with one of the early adaptations for the multi-objective scenario being proposed by Xue et al. (2004). As multi-objectives are considered, a dominance-based selection operator is used instead, the same as described in Definition 3.
In addition to a good approximation of the Pareto front, it is also often desirable for well-spread solutions throughout the objective-space, and not too clustered towards certain points. In the last decades, algorithms employed a variety of methods to achieve diversity and uniformity, with examples including function sharing (Srinivas & Deb, 1994), crowding , and ε-dominance (Laumanns et al., 2002).

Multi-objective differential evolution with spherical pruning (sp-MODE)
The sp-MODE, a variant of the multi-objective DE proposed in (Reynoso-Meza et al., 2010) and complemented in (Reynoso-Meza et al., 2014), is one such algorithm that focuses on the diversity of Pareto front solutions, working on fixing shortcomings of previous diversity-ensuring methods. The method employed by the algorithm, denominated Spherical Pruning, is based on allocating solutions on sectors of a hyper-sphere. For problems with three objectives, such as the one tackled in this work, the process can be visualized first as a projection of the normalized Pareto front approximation onto the first octant of the unit sphere. Considering a system of spherical coordinates ( r, , , ) the sphere surface is then divided into sectors by considering ∈ equal divisions for and , in the interval [0, 2 [ (Fig. 1). By allocating a single solution per sector, whichever happens, to be closer to the origin (a lesser value for r ), the algorithm can ensure diversity and well-spread solutions on the objective space.
In addition to the spherical pruning to ensure diversity, the algorithm also employs a secondary method, that ensures pertinency. The acquisition of a set of solutions with different compromises is the goal of a MOP, however, not all possibilities may be pertinent to the designer, but only a subset of them. The sp-MODE employs a physical programming approach to assign meaningful pertinency intervals for each of the sought objectives, in a manner more transparent than other methods, such as weighted sum and goal attainment .
Although the physical programming guides the algorithm towards the area of interest of the objective space, and the spherical pruning ensures well-spread solutions around that area, the latter can still be improved. As shown in Fig. 1, the constructed sectors on the sphere surface are highly non-uniform, with a different concentration of sectors depending on the location on the sphere surface. That can not only introduce a bias

Uniform spherical pruning
The spherical pruning is implemented in a way that allows that mechanism that constructs the spherical sectors to be extracted and defined as a function (Eq. (11)).
where a represents the k − 1 angles that define the point in hyper-spherical coordinates, and ∈ its hyper-spherical sector, defined by k − 1 values. Since the point is projected onto the hyper-sphere, the r parameter is not utilized.
The default sector-constructing function simply outputs a division of the inputs, resulting in the non-uniform grids, a result that the method described below focuses on correcting.
The method in question, proposed in (Roşca, 2010), is a sphere division method that maps points from a square ( a, b ), to a circumference ( X, Y ), to one of the hemispheres of a sphere ( x, y, z ), as displayed in Fig. 2.
Since the method maintains the areas between transformations, a uniform grid can be defined on the square (Fig. 2), and the reverse transformations (Eq. (13) and Eq. (14)) applied to find out where the points on the sphere will fall on the square. Based on the strategy for surface division described in this section, this study has been modified the original version of the sp-MODE proposed by Reynoso-Meza et al. (2010) by including a uniform grid proposed by Roşca (2010). The proposed approach creates a grid of cells of uniform areas on the surface of the sphere, representing the objective space, and consequentially, a uniform and less biased spherical pruning when compared to the original strategy presented by Reynoso-Meza et al. (2010).

Design parameters
As presented before, this work was developed to solve a real problem in the additive manufacturing process of product design using 3D printing technology. The experimental study consists of investigating the product made by Project 3000 Professional 3D Printer (3D Systems, 2011), from the conversion of the file for the Standard Triangle Language, or Standard Tessellation Language (STL), which represents a native exchange format among Stereolithography CAD systems intending to reduce the geometry to a triangular set of surface facets defined by three vertices and a normal vector. This representation is presented in the software screen illustrated in Fig. 3a.
The cylindrical specimen (Fig. 3b) was developed based on international standard ISO 4287:1997 (ISO, 1997) by the CAD software CATIA V5®. The 2 k factorial design, considering three input variables (k = 3) tested at two levels of control with three replicates each (2 3 × 3 = 24). This experimental study had limited material resources and experiments were carried out with a limited number of replicates. Thus, the way of working for this study was drawn as Zhang et al. (2009) and Chen and Zhao (2016), who obtained significant results using a limited number of replicates. Zhang et al. (2009) carried out a fractional factorial design of experiments applied to synthesis of ZSM-5 zeolite using 16 experiments, two replicates for each of the eight experiments. Chen and Zhao (2016) performed a design of experiments for binder blasting additives manufacturing process to optimize the layer thickness, printing saturation, heater power ratio, and drying time using two replicates for each combination of factors, obtaining a reliable result. For this experiment, epoxy resin 500 g cartridges with internal chip were assumed. The three input variables: (i) the layer thickness (LT), which defines the density of the piece; (ii) the deposition rate (DR), which sets the injection quantity of droplets per second; and (iii) the print orientation (O), which defines the angle of inclination of the part concerning to the deposition of the layers (see Table 1). The experimental study considered the length and thickness of the piece as noise parameters, due to the statistical analysis for these parameters, there is no statistically significant difference. It has been used the same extrusion head for the printing paths in all experimental runs. Thus, the noise parameters did not affect the output variables. The build direction used for building the pieces was horizontal axis (x-axis) in order to minimize the material losses, energy cost, and cycle time.
The output variables of the manufacturing process were roughness, dimensional accuracy, and mechanical resistance. The roughness was evaluated by the irregularities and protrusions above the surface. The unit of measure for roughness used for the specimens was the average roughness (Ra) given by the Mitutoyo SJ 201P model (Mitutoyo, 2004). When considering the dimensional accuracy, it was evaluated by the difference between the measured value and the reference value (Ø = 10 mm). The measurements were performed by instrument Micrometer Mitutoyo 0-25 mm range and resolution 0.01 mm (Mitutoyo, 2016). Regarding mechanical resistance, tensile testing was considered where test pieces were gripped in jaws and stretched by moving the grips apart at a constants rate while measuring the load (see Fig. 3). Measurements were conducted in a Lloyd dual-column tensile testing machine LD Series (AZO Materials, 2017), and the unit for that variable is N/mm 2 .

Comparison of metrics
To measure how well the modified algorithm will perform over the proposed problem, three Pareto comparison metrics are utilized to evaluate both the convergence and diversity of the obtained Pareto approximations for each algorithm, the: (i) hyper-volume, (ii) inverted generational distance, and (iii) g-Indicator.
The hyper-volume (H), proposed by Zitzler and Thiele (1999), is a metric focused on evaluating the convergence of the approximated Pareto front. It is defined as the dominated area "behind" the Pareto front, up until a point of reference (usually the NADIR). The hyper-volume is a unary metric that generates a larger value the closer the approximated front is to the true Pareto (Fig. 4a). Table 1 Variables statement for experiments *Once accuracy objective fluctuates around "0" with positive and negative values, a generalization must be considered to turn it a minimization problem, then it has been adopted Accuracy =|Accuracy| The g-Indicator (G), which was proposed by Lizárraga et al. (2008), is a metric to evaluate the diversity and distribution of solutions on the approximated front. In its generalized form, the metric is defined by the union of the hyper-volume covered by hyper-spheres of radius U centered at each point of the Pareto front approximation. The g-Indicator is a m metric that calculates its value for m fronts simultaneously, with the radius U being calculated as half the average distance of all solutions, of all fronts, to their closest neighbor (Eq. (15)). Additionally, the process is simplified by projecting the given approximated fronts to a lower dimension, such as a three-dimensional front to a two-dimensional plane (Fig. 4b).
where m represents the number of approximated fronts, and r ij the distance of the i-th solution of the j-th front, to its closest solution.
The inverted generational distance (IGD), proposed by Coello Coello and Sierra (2004), is a modified version of the classical generational distance (GD) used to evaluate both convergence and diversity evaluation of a given approximated front, in comparison to a reference front. The original GD uses an average of the distances of every solution on the approximated front, to the closest solution on the reference, accounting solely for convergence. The IGD uses the distances from the solutions on the reference to the approximated front, considering how uniformly they are distributed. As the true Pareto is usually unknown for a real-world problem, it has been applied the aggregated approximated fronts from all runs as reference.

Analysis of experiments
This section describes the statistical analysis of the experiments presented in the previous section. Firstly, a statistical ANOVA analysis was conducted in a way to establish the problem over the 3D printing technology process. After concluding the prototype, it was possible to establish the contribution of the input variables in each one of the outputs. Therefore, the outputs (roughness, accuracy, and resistance) were analyzed regarding the input variables (layer thickness, deposition rate, and print orientation). Starting with the analysis of roughness, the average roughness found in this case through the software was 43.83 Ra. Table 2 shows the analysis of variance for roughness. All main effects and two-way interactions (thickness × rate and thickness × orientation) are statistically significant (p < 0.05), indicating the difference between levels of the input variables with influence in the roughness output. When the p-value for the statistical F-test of the overall significance test is less than your significance level, it is possible to reject the null hypothesis and conclude that the model provides a better fit than the intercept-only model (Snedecor & Cochran, 1989).   Figure 5a presents the percentage of contribution of every input variable in the roughness. By using a Half-Normal plot of factor effects highlighting that, the main effects are provided by layer thickness, deposition rate, orientation, and the two-way interactions between them. The Half-Normal plot of factor effects allows one to detect the factor and interaction that are most important to the process, where an optimization design could be applied. This plot shows the absolute values of the effects and draws a reference line on the chart (dashed red line). Any effect that goes beyond this line is considered potentially important (Sant'Anna, 2015). Figure 5b shows the cube plot with adjusted averages for roughness. Cube plots display the average response values at all combinations of process and design parameter settings. In this way, it is possible to determine the best and the worst combinations of factor levels with the purpose to achieve the desired optimum response (Myers et al., 2016). In the sequence, it was possible to define the equation that conducts to roughness values in terms of the set of inputs, being so:   In the same way, as performed for roughness, the dimensional accuracy variable was evaluated, and the average was 0.0387 mm. Table 3 shows the analysis of variance considering the same approach that was adopted for roughness. Figure 6a shows the Half-Normal plot with the percentage of the contribution that every input variable has in the dimensional accuracy, highlighting that the main effects are caused by layer thickness, deposition rate, the two-way interactions between them, and the two-way interactions between layer thickness and orientation. Figure 6b shows the cube plot with adjusted averages for dimensional accuracy, highlighting the effects provided by layer thickness, deposition rate, and orientation. By applying the regression technique, it was possible to define the equation that conducts to dimensional accuracy values in terms of inputs:

Roughness
The mechanical resistance variable was also analyzed with an average of 260.6 N/mm 2 . Table 4 shows the analysis of variance for mechanical resistance, while Fig. 7a reports a Half-Normal plot with the percentage of the contribution that each input has in mechanical resistance. In this case, the orientation became the most relevant variable in this   Figure 7b shows the cube plot with adjusted averages for mechanical resistance. By assuming the same regression technique previously mentioned, the equation that defines the mechanical resistance in terms of the inputs can be obtained: Statistical analysis allows obtaining the optimal setting considering all the input parameters. By assuming the analysis of variance test, half-normal, and cube plots were possible to obtain a solution for this problem. However, this approach works well in problems that involve a low quantity of variables in decision-making. When the complexity of the problem significantly increases relate to the precision of parameters, as presented in the current problem described in this work, where it is possible to combine three possible inputs with three affected outputs, the solution in most cases is not trivial. Figure 8a-c shows the surfaces responses plot for Roughness, Accuracy, and Resistance models based on the inputs. Note that for orientation of printing, it was only considered the at 0° grad once the results reached for 90° grad of printing orientation were quite poor. This decision also supports having a 3D visualization on surfaces responses, otherwise, it should be a bi-cubic graph, which is non-trivial to visual analysis.

Algorithm evaluations
Based on the previous assumption, the usp-MODE algorithm has been applied in the case study to illustrate its application to more complex cases. In this way, the usp-MODE supports both search process and decision-making to find the best solution based on the compromise between inputs and outputs, as infinite possibilities can be assumed for the design variables within the predefined search space, where {x 1 , x 2 , x 3 }. Additionally, to ensure a bias-free evaluation of the proposed algorithm, it is applied to multiple benchmark functions, as well as the DOE problem described earlier. This is done to ensure an acceptable performance under Pareto sets and fronts of different shapes and characteristics. The benchmark functions are presented in Cheng et al. (2018), which describes a methodology for the evaluation of algorithms for many-objective optimization problems. These functions, although meant to be used in problems with four objectives or more, are defined for three-objective problems.
By assuming 24 real experiments to set the equation that establishes the relation between inputs and outputs, the intention to apply the usp-MODE algorithm is to expand the search space of the input variables that up to now were considered just in two levels each. The purpose of the problem is through the three inputs (LT, DR, and O) to minimize both roughness and accuracy, at the same time maximizing the mechanical resistance. Finally, considering the Roughness, Accuracy, and Resistance models (by Eqs. (9-11)) for applying the usp-MODE algorithm, it was possible to obtain values for inputs to obtain a quasi-optimal configuration for the outputs.
Considering that inputs belong to the machining process, it was decided to keep the orientation only to 0º or 90º values given the difficulty to use different angles. On the other hand, both layer thickness and deposition rate were searched inside a range with boundaries determined by the values shown in Table 5, being so: LT in the interval [0.032, 0.064] µm and DR between [150, 300] u/sec have been defined.
The parameterization of the NSGA-II was performed base on values found in the literature (Alvarado-Iniesta et al. 2017), Table 5 shows the complete list of parameters applied. For a fair comparison, both the population and the iterations of the remaining two algorithms (sp-MODE and usp-MODE) were chosen to result in the same number of function evaluations. The remaining parameters were kept as in the default implementation suggested on the original paper, including the number of arc divisions equal to 10 times the number of optimization objectives.
To account for the stochastic nature of the algorithms, 31 individual runs were performed. The three metrics were calculated for each of the benchmark runs, and their  Table 9. Firstly, the results of the H metric showed an improvement in convergence regarding the proposed modification, albeit small. As the variations were not large, it could be verified that the modification performed in the technique to improve diversity did not impact the algorithm's convergence.
The results of the G metric displayed more considerable improvements (function 1, 6, 8 and 13), if compared to the previous case of the H metric. In the same fashion, for the cases where the other two algorithms outperformed the usp-MODE, the improvements were considerably less significant. Additionally, results for the proposed algorithm showed a smaller standard deviation, in most cases. The IGD metric accounts for the number of solutions by using an average for the distances. Although not as significant as the previous metric, the usp-MODE still showed improvements over the default sp-MODE. In the same way as the last metric, improvements were considerable, apart from the 6th benchmark function, due to the degenerate property of the Pareto front.
Results for the optimization itself show similar behavior with the previous three metric benchmarks. While overall improvements were observed in all three metrics, the hypervolume displayed the less significant one (5.73% over NSGA-II, and 1.23% over sp-MODE). The results for the G metric showed the most significant improvement (194.34% over NSGA-II, and 32.20% over sp-MODE), but also a higher variance. Finally, although the results for the IGD metric were less pronounced (63.38% over NSGA-II, and 20.59 over sp-MODE), they also showed the smallest variance. These results not only show that the convergence capacity of the algorithm was not compromised by the change on the spherical pruning mechanism but the capability of reaching a more diversified Pareto front was achieved.
The aggregated approximated Pareto front of all runs for the three algorithms was constructed, for better visualization of different compromises between objectives, and it is presented in Fig. 9.

Discussions
The proposed approach used to multiobjective optimization problem in this real 3D printing process was based on the integration of the design of experiments with optimization algorithm, optimizing three conflicting objectives in a suitable finite-size illustrated by Pareto front, once it offers a decision-making tool with a good result of the optimal design. Myers et al., (2016) recommend performing a statistical test to define if there is systematic nonlinearity curvature present. Due to that, authors have performed the statistical test to analyze the lack of fit to first-order model with interaction showed by Eqs. (16-18) from the designed experiment. The linearity was statistically significant indicating that the straight-line fit was very satisfactory as can be seen in the half-normal probability plot of the effects by Figs. 5a, 6a, and 7a for the roughness, accuracy, and resistance variables. The overall proposed usp-MODE optimization algorithm of the parameters in 3D printing technology was achieved based on the following optimum values to input variables: layer thickness at low level (0.032 µm), deposition rate at the high level (300 u/sec), and print orientation at low level (0° grad). Another important result of this work was to verify the optimal values to output variables: surface roughness 33.33 Ra, dimensional accuracy 0,12 mm, and mechanical resistance 204.40 N/mm 2 , improving the outputs performances of the machining process. This study highlights the relevance of an experimental optimization for the machining parameters process based on planned experimental tests and a robust optimization algorithm, that allows to minimize the loss of material and operating cost.
The proposed optimization approach's results are corroborated with the literature, Bhavsar et al. (2015) assert that the NSGA-II algorithm optimization generated optimal solutions of Pareto front for conflicting output variables performance, as material removal rate and surface roughness in FIB micro-milling process. Rao et al. (2017) showed that multiobjective optimization proved optimal solutions for multiple conflicting objectives as cutting velocity, surface quality, and dimensional accuracy to improve the performance of the machining processes. Alvarado-Iniesta et al. (2019) assert that the NSGA-III algorithm presented reliable results of Pareto fronts of multiobjective optimization problem used in the plastic injection molding process.

Conclusions
This research proposes a hybrid approach to optimize the parameters, integrating both design of experiment and multiobjective optimization methods. Based on the proposed hybrid approach, a support solution for 3D printing technology as an alternative to classical parametrization design used in the machining process, where the proposed algorithm usp-MODE was considered to select the "optimal" set of parameters in 3D printing technology. By assuming three input variables (layer thickness, deposition rate, and print orientation), and three outputs (surface roughness, dimensional accuracy, and mechanical resistance), the influence of each input variable, and the combination of them, to improve the process outputs performances were evaluated.
This study applied the 2 k factorial design, which provided reasonable results as important information about the influence of both individual and interaction variables in each one of the outputs. In this way, it was possible to define the best set for the input variables that fit the machining process, can be applied to distinct 3D printing technology processes. This procedure permits the extension of the DOE method adopted in this study to more complex cases where a considerable number of variables are considered. In this case, it is possible to increase the complexity of the problem without changing the approach presented in this study.
In the search procedure, it has been proposed the usp-MODE approach, and a benchmark case with NSGA-II and the classical sp-MODE algorithm has been carried out to show the benefits and gains of the proposed algorithm. It can be seen the proposed usp-MODE is a promising technique for solving many machining process optimization problems, providing a wider and more diverse set of Pareto solutions in most cases.
Future works will focus on the evaluation of distinct evolutionary algorithms applied to the multiobjective optimization problems reported in this study, to compare results with the use-MODE approach. Authors have also considered the possibility of applying the usp-MODE algorithm to different experiments meaning many-objectives problems, in this case, the reliability and robustness of the proposed algorithm might be evaluated in different ways and approaches. Finally, to bring a full analysis of the experimental design approach, other inputs could be considered, as energy efficiency and costs.