A Consistent Procedure Using Response Surface Methodology to Identify Stiffness Properties of Connections in Machine Tools

Accurate finite element models of mechanical systems are fundamental resources to perform structural analyses at the design stage. However, uncertainties in material properties, boundary conditions, or connections give rise to discrepancies between the real and predicted dynamic characteristics. Therefore, it is necessary to improve these models in order to achieve a better fit. This paper presents a systematic three-step procedure to update the finite element (FE) models of machine tools with numerous uncertainties in connections, which integrates statistical, numerical, and experimental techniques. The first step is the gradual application of fractional factorial designs, followed by an analysis of the variance to determine the significant variables that affect each dynamic response. Then, quadratic response surface meta-models, including only significant terms, which relate the design parameters to the modal responses are obtained. Finally, the values of the updated design variables are identified using the previous regression equations and experimental modal data. This work demonstrates that the integrated procedure gives rise to FE models whose dynamic responses closely agree with the experimental measurements, despite the large number of uncertainties, and at an acceptable computational cost.


Introduction
Machine tools are stationary, power-driven industrial devices used to manufacture workpieces under user and technological requirements. The most demanded requirements are accuracy and precision, which mainly depend on the static deformation and dynamic behavior of the machine tool under variable cutting forces. Assembly errors, tool trajectory errors, and the effect of thermal sources are also important issues [1]. Therefore, machine tool manufacturers devote strong efforts to perform the appropriate static, modal, and dynamic analyses of the machines, in order to determine the stresses and displacements, natural frequencies, and mode shapes. The final aim is to identify and analyze the vibration sources under different operating conditions, in order to minimize their effects on the surface finishing of the workpieces, stop the appearance of regenerative vibrations or chatter, and slow down the swift wear of the tools [2,3].
Today, the design process of modern machine tools is developed under virtual environments, where the finite element method (FEM) is widely used and particularly advised. The FEM provides a discretized model of the machine tool, whose purpose is to reproduce the real behavior of the structure. linear axes and is made up of four main modules, namely, a bed frame, column, framework, and ram, which slide over roller type linear guideways. Two servo motors, directly coupled to ball-screws supported by bearings at both ends, provide the displacement along Y-and Z-axes, while the movement in the X-axis is performed by a linear motor. The machine is joined to a concrete basement by anchor bolts and leveling elements adjust and align the bed frame.
Firstly, a FE model of the machine tool has been defined ( Figure 1), which consists of 12,804 nodes and 14,983 elements, mainly shell and solid brick elements. The connections between the different components and the connection to the foundation have been modeled by linear spring elements. In this way, the contact elements and friction coefficients in the FE formulation are avoided, reducing the complexity and keeping the model linear. These linear springs characterize the previously mentioned linear guideways, ball-screws, and bolts, and are incorporated into the FE model in their locations. The anchor bolts connecting the bed frame to the basement behave rigidly, so high stiffness spring elements have been used in the modeling process. Linear guideways have been modeled assigning average stiffness values in two directions, perpendicular and transverse to the direction of movement, based on the stiffness curves provided by the guideway supplier [27], and very low stiffness values along the directions where the movement is developed. A similar modeling has been followed for the ball-screws, although, in this case, low stiffness values have been set in perpendicular and transverse directions to the main movement [2,28,29]. Moreover, the tool holder has been simulated as a beam, and spindle, servo motors and the face milling cutter as lumped masses. Finally, solid brick elements have been used to model the primary and secondary sections of the linear motor with a spring element between them. directly coupled to ball-screws supported by bearings at both ends, provide the displacement along Y-and Z-axes, while the movement in the X-axis is performed by a linear motor. The machine is joined to a concrete basement by anchor bolts and leveling elements adjust and align the bed frame.
Firstly, a FE model of the machine tool has been defined ( Figure 1), which consists of 12,804 nodes and 14,983 elements, mainly shell and solid brick elements. The connections between the different components and the connection to the foundation have been modeled by linear spring elements. In this way, the contact elements and friction coefficients in the FE formulation are avoided, reducing the complexity and keeping the model linear. These linear springs characterize the previously mentioned linear guideways, ball-screws, and bolts, and are incorporated into the FE model in their locations. The anchor bolts connecting the bed frame to the basement behave rigidly, so high stiffness spring elements have been used in the modeling process. Linear guideways have been modeled assigning average stiffness values in two directions, perpendicular and transverse to the direction of movement, based on the stiffness curves provided by the guideway supplier [27], and very low stiffness values along the directions where the movement is developed. A similar modeling has been followed for the ball-screws, although, in this case, low stiffness values have been set in perpendicular and transverse directions to the main movement [2,28,29]. Moreover, the tool holder has been simulated as a beam, and spindle, servo motors and the face milling cutter as lumped masses. Finally, solid brick elements have been used to model the primary and secondary sections of the linear motor with a spring element between them. Table 1 describes the main parameter values of the FE model.    Then, using the Lanczos solver, the free motion of the structure has been analyzed by calculating the natural frequencies and mode shapes from the assembled mass and stiffness matrices of the numerical model. As the connection between the bed frame and the basement is considered in the FE model, these modal parameters correspond to the in situ configuration of the machine. According to several tests developed under chatter conditions [30], the frequency range of interest has been defined as 10 Hz to 120 Hz. The natural frequencies are shown in Table 2.

Experimental Modal Analysis
In order to experimentally determine the dynamic characteristics of the machining center, a multiple reference impact test was performed, using, as references, point 5 along X and Y directions and by exciting the system with an instrumented hammer. The translational acceleration responses in the X-, Y-, and Z-axes were measured in 75 points using triaxial accelerometers, so the accelerance frequency response functions (FRFs) corresponding to 225 degrees of freedom were obtained. The total number of measured FRFs was 450. Figure 2 shows a wire frame model representation of the test structure. The references are identified with arrows. Then, using the Lanczos solver, the free motion of the structure has been analyzed by calculating the natural frequencies and mode shapes from the assembled mass and stiffness matrices of the numerical model. As the connection between the bed frame and the basement is considered in the FE model, these modal parameters correspond to the in situ configuration of the machine. According to several tests developed under chatter conditions [30], the frequency range of interest has been defined as 10 Hz to 120 Hz. The natural frequencies are shown in Table 2.

Experimental Modal Analysis
In order to experimentally determine the dynamic characteristics of the machining center, a multiple reference impact test was performed, using, as references, point 5 along X and Y directions and by exciting the system with an instrumented hammer. The translational acceleration responses in the X-, Y-, and Z-axes were measured in 75 points using triaxial accelerometers, so the accelerance frequency response functions (FRFs) corresponding to 225 degrees of freedom were obtained. The total number of measured FRFs was 450. Figure 2 shows a wire frame model representation of the test structure. The references are identified with arrows. From the measured FRFs, a polyreference Least Squares Complex Frequency (pLSCF) estimator was used to extract the system modal parameters. Table 3 shows the natural frequencies and a brief description of the different mode shapes. From the measured FRFs, a polyreference Least Squares Complex Frequency (pLSCF) estimator was used to extract the system modal parameters. Table 3 shows the natural frequencies and a brief description of the different mode shapes.

Comparison between FE and Experimental Modal Data
At this point, there are two sets of different results, related to the numerical and experimental models. The next step is to evaluate the correspondence between them, because it is necessary that both models show a considerable degree of correlation, in order to improve the FE model successfully.
Firstly, the geometrical correlation has been developed to match the different coordinate and unit systems used in the models, and, then, the mode shape correlation has been performed to establish a reliable pairing between the numerical and experimental modes. A very useful indicator to compare and contrast the modal vectors from the different sources is the modal assurance criterion (MAC) [31]. The MAC shows the degree of linearity between two modal vectors, ϕ FEA (FEA-finite element analysis) and ϕ exp , as follows: The MAC can take on values from 0, showing a lack of correspondence between the modal vectors, to 1, which means that the modal vectors are the same. Table 4 shows the frequency differences and MAC values between the FE and experimental responses, wherein the MAC values corresponding to the paired mode shapes have been bolded. These correlation results can be considered sufficient for a large number of practical applications [32,33], as the mean frequency difference is 3.3%, and the mean MAC value is 87.0%. Nevertheless, there are still moderate differences between several natural frequencies, which confirm that it is necessary to improve the FE model. Therefore, the main goal of the following procedure will be to match the numerical frequencies to the experimental ones, while maintaining or improving the MAC pairing values.

Two-Level Designs for Parameter Screening
Among the different types of experimental designs [14], the factorial designs are widely used to identify, at the initial stages, the main variables that affect any process or system (i.e., as screening experiments). The basic design is a two-level or 2 k design, where k is the number of variables and each of them takes an upper and a lower level. A complete trial of such a design needs 2 k runs and allows for estimating the linear effects of the k variables and their interactions.
Nevertheless, as the number of variables, k, increases, the number of runs in the trial also increases, but dramatically, and interactions between three, four, and more variables appear. Assuming that the highest interactions are negligible, it would be possible to obtain information concerning the effects of the variables and low-order interactions by running a part or fraction of the complete factorial design, 2 k−p , where p indicates the fraction chosen (1/2 p ). The so-called resolution V design is especially interesting, which provides information about the contribution of variables and two-factor interactions, mixed with higher-order interactions. As these are negligible, the fractional designs are better than the complete factorial designs, because the number of runs diminishes considerably.
Once the trial has been finished, the next step is to identify the significant factors and interactions by performing an analysis of variance on the results. According to ANOVA, the variability of the results in an experiment that is dependent on several variables, is the sum of variability due to each factor, plus that contributed by the interaction between the factors, and that added by the internal error. Also, using ANOVA, the sum of squares (SS) can be used as a measure of the overall variability, so that the greater the SS due to a factor, the larger its importance on the process or system. Thus, it will be possible to find out which variables and interactions are the most significant.

Response Surface Methodology to Develop an Optimal Mathematical Model
The purpose of response surface methodology is to build an explicit function to approximate the actual relationship between the variables, x i , and a response, y, involved in an engineering problem. That function, preferably a low-order polynomial, is in fact a regression model, less expensive to evaluate, which can be used to predict the response developed in the system under a specific combination of variables.
In general, the behavior of the industrial processes and mechanical systems cannot be explained by linear functions [19,23,34], so, in the following section, the second-order models (Equation (2)) and the experimental designs that are preferable to adequately estimate these models will be examined.
In Equation (2), β 0 is the average value of response; y, β i , β ii , and β ij are the partial regression coefficients; ε is the error term; and k is the number of variables.
One of the most popular designs for fitting second-order models is the central composite design (CCD). It is built in a sequential way, based on a two-level factorial (2 k ) design, plus 2k axial and n C center points. The points added to factorial design allow an efficient estimation of the possible curvature of the model.
Firstly, a set of responses y is obtained on the completion the experiments of the central composite design. Then, the values of these responses and design variables are substituted in Equation (2), and rewritten in matrix form as follows: Equation (3) is solved using the least squares method, by minimizing the sum of the squares of the errors ε i . That leads to a least squares estimator of β, as follows: At this point, an initial second-order model is completely defined using all of the design variables and interactions. Then, it is necessary to perform an analysis of variance to check the significance of each parameter and the adequacy of the regression model. For this last purpose, various statistical parameters can be used, such as, the coefficient of determination R 2 , the adjusted R 2 , and the predicted R 2 [17]. These coefficients are all expected to be close to 1.0, which would mean that the regression model, y RSM , explains the response, y, properly and that it also predicts adequately new responses.
Nevertheless, if there are substantial differences between them, the least significant parameter is removed using the t-test, and a new regression model, Equation (2), is built and the analysis is repeated until the remaining parameters are all significant. On the completion of the iteration process, the optimum response surface model can be considered as adequate to carry on the next stage of the improvement procedure.

Identification of Updated Values of the Design Variables using the Optimum RS Model
Once the mathematical relationships between the design variables and responses have been established, the final step is to identify those values of the design variables that lead to the responses that better fit the experimental ones. This is actually an inverse multi-objective constrained optimization problem, and nonlinear programming techniques can be used to solve it.
Another alternative approach is based on the so-called desirability function [35], which is explained in the following. Firstly, each estimated response, y RSMi , is turned into a desirability function, d i , as follows: d i = 0, y RSMi < y LOWi and y UPi < y RSMi (7) where y OBJi is the target experimental response, and y LOWi and y UPi are the lower and upper limits for each response ( Figure 3).
Materials 2018, 11, x FOR PEER REVIEW 8 of 21 b = (X T · X) −1 · X T · y (4) At this point, an initial second-order model is completely defined using all of the design variables and interactions. Then, it is necessary to perform an analysis of variance to check the significance of each parameter and the adequacy of the regression model. For this last purpose, various statistical parameters can be used, such as, the coefficient of determination R 2 , the adjusted R 2 , and the predicted R 2 [17]. These coefficients are all expected to be close to 1.0, which would mean that the regression model, yRSM, explains the response, y, properly and that it also predicts adequately new responses.
Nevertheless, if there are substantial differences between them, the least significant parameter is removed using the t-test, and a new regression model, Equation (2), is built and the analysis is repeated until the remaining parameters are all significant. On the completion of the iteration process, the optimum response surface model can be considered as adequate to carry on the next stage of the improvement procedure.

Identification of Updated Values of the Design Variables using the Optimum RS Model
Once the mathematical relationships between the design variables and responses have been established, the final step is to identify those values of the design variables that lead to the responses that better fit the experimental ones. This is actually an inverse multi-objective constrained optimization problem, and nonlinear programming techniques can be used to solve it.
Another alternative approach is based on the so-called desirability function [35], which is explained in the following. Firstly, each estimated response, yRSMi, is turned into a desirability function, di, as follows: where yOBJi is the target experimental response, and yLOWi and yUPi are the lower and upper limits for each response ( Figure 3). Then, a global desirability function, D, is built as the geometric mean of individual desirabilities, di, as follows: where u is the total number of the experimental responses. Then, a global desirability function, D, is built as the geometric mean of individual desirabilities, d i , as follows: where u is the total number of the experimental responses. Finally, the results are ranked in decreasing desirability order and the values of the design variables that maximize the global desirability D are selected.

Initial Selection of Candidate Design Variables
In order to improve the FE model, first, it is necessary to select the design variables to work with. There are a large number of design parameters to be considered in this machining center, but, in fact, the main uncertainties in the FE model are concentrated on connections, as follows: Stiffness values of the connection elements between main components of the machine tool ( Figure 4); 2.
Stiffness values assigned to the elements that attach the machine tool to the foundation ( Figure 5); and 3.
Stiffness value along X direction of the connection element, between the primary and secondary sections of the linear motor.
Materials 2018, 11, x FOR PEER REVIEW 9 of 21 Finally, the results are ranked in decreasing desirability order and the values of the design variables that maximize the global desirability D are selected.

Initial Selection of Candidate Design Variables
In order to improve the FE model, first, it is necessary to select the design variables to work with. There are a large number of design parameters to be considered in this machining center, but, in fact, the main uncertainties in the FE model are concentrated on connections, as follows:  Nevertheless, the number of variables is still large, three stiffness values for the joints to the foundation, the stiffness value for the inner connection of the linear motor, and six stiffness values for the connections between the modules of the machining center (Table 5). Therefore, first, it is necessary to determine which variables affect, to a large extent, each model response and, hence, remove those whose influence is negligible. For this purpose, a resolution V design would be the most convenient (i.e., in our case, a 2 10−3 fractional design with 128 runs) [14]. However, it is still too laborious to manage such a number of runs. Furthermore, some parameter combinations could lead to inappropriate model responses, due to the presence of the constraints among them. Thus, in order to facilitate the analysis and gain a progressive comprehension of the significance of each design variable and interaction, instead of a single 2 10−3 fractional design, an alternative trial with seven 2 5−1 designs (16 runs each) has been performed ( Table 5). Each variable has been paired up with the rest Finally, the results are ranked in decreasing desirability order and the values of the design variables that maximize the global desirability D are selected.

Initial Selection of Candidate Design Variables
In order to improve the FE model, first, it is necessary to select the design variables to work with. There are a large number of design parameters to be considered in this machining center, but, in fact, the main uncertainties in the FE model are concentrated on connections, as follows:  Nevertheless, the number of variables is still large, three stiffness values for the joints to the foundation, the stiffness value for the inner connection of the linear motor, and six stiffness values for the connections between the modules of the machining center (Table 5). Therefore, first, it is necessary to determine which variables affect, to a large extent, each model response and, hence, remove those whose influence is negligible. For this purpose, a resolution V design would be the most convenient (i.e., in our case, a 2 10−3 fractional design with 128 runs) [14]. However, it is still too laborious to manage such a number of runs. Furthermore, some parameter combinations could lead to inappropriate model responses, due to the presence of the constraints among them. Thus, in order to facilitate the analysis and gain a progressive comprehension of the significance of each design variable and interaction, instead of a single 2 10−3 fractional design, an alternative trial with seven 2 5−1 Nevertheless, the number of variables is still large, three stiffness values for the joints to the foundation, the stiffness value for the inner connection of the linear motor, and six stiffness values for the connections between the modules of the machining center (Table 5). Therefore, first, it is necessary to determine which variables affect, to a large extent, each model response and, hence, remove those whose influence is negligible. For this purpose, a resolution V design would be the most convenient (i.e., in our case, a 2 10−3 fractional design with 128 runs) [14]. However, it is still too laborious to manage such a number of runs. Furthermore, some parameter combinations could lead to inappropriate model responses, due to the presence of the constraints among them. Thus, in order to facilitate the analysis and gain a progressive comprehension of the significance of each design variable and interaction, instead of a single 2 10−3 fractional design, an alternative trial with seven 2 5−1 designs (16 runs each) has been performed ( Table 5). Each variable has been paired up with the rest at least once, so that after the completion of the whole set of 2 5−1 experiments, it has been possible to look into the effects of all of the design variables and two-factor interactions by means of ANOVA. Prior to conducting the fractional designs, the range of each variable was decided (Table 5) according to load-deformation curves [27] and previous works [30].
On completion of the trial, the total corrected sum of squares, SS T (Figure 6), and the sum of squares of each factor and the two-factor interactions, mixed with higher interactions (SS i and SS ij ) (Figures 7 and 8), have been obtained as a measure of the variability, for all of the frequencies and MAC values. Firstly, for each response, the SS T has been examined, as some designs add much more variability than others, because of the variables involved. Thus, Figure 6 shows that MAC 1 and MAC 2 are not affected by the changes in the design variables, and that the variability of f FEA2 is negligible. As this frequency matches its experimental pair (Table 4), it has been omitted in later analyses. at least once, so that after the completion of the whole set of 2 5−1 experiments, it has been possible to look into the effects of all of the design variables and two-factor interactions by means of ANOVA. Prior to conducting the fractional designs, the range of each variable was decided (Table 5) according to load-deformation curves [27] and previous works [30]. On completion of the trial, the total corrected sum of squares, SST (Figure 6), and the sum of squares of each factor and the two-factor interactions, mixed with higher interactions (SSi and SSij) (Figures 7 and 8), have been obtained as a measure of the variability, for all of the frequencies and MAC values. Firstly, for each response, the SST has been examined, as some designs add much more variability than others, because of the variables involved. Thus, Figure 6 shows that MAC1 and MAC2 are not affected by the changes in the design variables, and that the variability of fFEA2 is negligible. As this frequency matches its experimental pair (Table 4), it has been omitted in later analyses. In order to determine which variables and two-factor interactions provide the largest influence on the variability of natural frequencies and mode shapes, SSi and SSij have been gradually analyzed in each 2 5−1 design. As an example, Figures 7 and 8 illustrate the percentage contribution of each parameter in design number 4.   kX210, kY3, kZ4, kX8, and kY9, respectively; 6-15, two factor interactions kX21-kY3, kX210-kZ4, kX210-kX8, kX210-kY9, kY3-kZ4, kY3-kX8, kY3-kY9, kZ4-kX8, kZ4-kY9, and kX8-kY9, respectively.  kX210, kY3, kZ4, kX8, and kY9, respectively; 6-15, two factor interactions kX210-kY3, kX210-kZ4, kX210-kX8, kX210-kY9, kY3-kZ4, kY3-kX8, kY3-kY9, kZ4-kX8, kZ4-kY9, and kX8-kY9, respectively. From Figure 7, it is found that the design variables (1-5) have a larger influence on the natural frequencies than the two-factor interactions (6)(7)(8)(9)(10)(11)(12)(13)(14)(15). For instance, design variable 2 (kY3) dominates the 1st and 2nd natural frequencies, design variable 4 (kX8) has a huge influence on the 6th natural frequency and an important weight on the 3rd and 4th natural frequencies, and design variable 5 (kY9) governs the 5th natural frequency. However, only parameter 8 (interaction kX210-kX8) slightly affects the 3rd and 4th natural frequencies. Also, design variable 3 (kZ4) does not seem to have a notable influence on any natural frequency.
On the other hand, from Figure 8, some of the interactions play significant roles in MAC values (3, 4, and 6, mainly). In fact, MAC4 is heavily affected by parameter 8 (interaction kX210-kX8), which also provides an important contribution to the variance of MAC3. Also, parameter 11 (interaction kY3-kX8) causes approximately 35% of the total variability to MAC6. Nevertheless, in general, the individual design variables have a greater influence than the interactions. This analysis has been repeated for each 2 5-1 design, so that it has been possible to gradually gain a better insight into the influence of the design variables and two-factor interactions on the responses. Finally, those factors providing more than 99% of the total variability for the frequencies, and 97.5% for the MAC values have been selected (Tables 6 and 7) to continue with the improvement procedure.
In order to determine which variables and two-factor interactions provide the largest influence on the variability of natural frequencies and mode shapes, SS i and SS ij have been gradually analyzed in each 2 5−1 design. As an example, Figures 7 and 8 illustrate the percentage contribution of each parameter in design number 4. From Figure 7, it is found that the design variables (1-5) have a larger influence on the natural frequencies than the two-factor interactions (6)(7)(8)(9)(10)(11)(12)(13)(14)(15). For instance, design variable 2 (k Y3 ) dominates the 1st and 2nd natural frequencies, design variable 4 (k X8 ) has a huge influence on the 6th natural frequency and an important weight on the 3rd and 4th natural frequencies, and design variable 5 (k Y9 ) governs the 5th natural frequency. However, only parameter 8 (interaction k X210 -k X8 ) slightly affects the 3rd and 4th natural frequencies. Also, design variable 3 (k Z4 ) does not seem to have a notable influence on any natural frequency.
On the other hand, from Figure 8, some of the interactions play significant roles in MAC values (3, 4, and 6, mainly). In fact, MAC 4 is heavily affected by parameter 8 (interaction k X210 -k X8 ), which also provides an important contribution to the variance of MAC 3 . Also, parameter 11 (interaction k Y3 -k X8 ) causes approximately 35% of the total variability to MAC 6 . Nevertheless, in general, the individual design variables have a greater influence than the interactions. This analysis has been repeated for each 2 5-1 design, so that it has been possible to gradually gain a better insight into the influence of the design variables and two-factor interactions on the responses. Finally, those factors providing more than 99% of the total variability for the frequencies, and 97.5% for the MAC values have been selected (Tables 6 and 7) to continue with the improvement procedure. Table 6. Summary of variables affecting responses.

Connection
Design Table 7. Summary of two-factor interactions affecting responses.

Responses Interactions
Several conclusions can be inferred from Tables 6 and 7, as follows: • The fractional factorial experiments have allowed for finding out the design variables and interactions that affect the responses. Therefore, the screening experiment has satisfactorily achieved the initial goal. • Two design variables do not influence the natural frequencies, namely the stiffness k X11 and horizontal stiffness k X21 of the connections to the foundation. • Three design variables are only significant for one natural frequency (f FEA5 ), namely the transverse stiffness k Z63 between the bed frame and foundations, and two stiffnesses between the modules of the machine tool, k Z13 and k Y9 . • MAC 5 and MAC 6 are affected by the largest number of interactions. Furthermore, some of them include design variables that do not influence them individually, for example, interaction k X11 -k X21 . This situation only appears in these two responses. In addition, the total number of design variables, considered both individually and in interactions, which affect each of these responses is nine (i.e., almost all). Nevertheless, along the complete set of fractional designs, the MAC 5 values were always larger than 80% and the MAC 6 values ranged from 68% to 73%. Thus, it has been decided not to carry on with the study of these responses, because the number of involved variables would lead to a costly analysis in the next step, while the benefits would be quite poor.

•
The natural frequencies f FEA3 and f FEA4 and the corresponding MAC values depend on the same group of design variables, k X210 , k Y3 , k Z4 , and k X8 . In addition, the natural frequency f FEA6 is dependent on three of these variables, k Y3 , k Z4 , and k X8 . Therefore, in the next step of the improvement process, these three natural frequencies will be analyzed together, so as to reduce the number of experiments necessary to define their meta-models. In addition, it is interesting to note that the mode shapes associated to these frequencies take place in plane XZ.

•
The natural frequency f FEA5 is affected by four variables that do not have any influence on the frequencies f FEA3 , f FEA4 , and f FEA6 , and, vice versa, the variables that affect these three frequencies do not provide any variability to the natural frequency f FEA5 . Moreover, some of the design variables representing stiffness in the X direction, k X8 and k X210 , do not affect the 1st and 5th mode shapes, whose principal movement is in plane YZ. Thus, it is concluded that the design variables are working collectively.

Development of Explicit Relationships between Design Variables and Responses
The next step of the improvement procedure is the definition of the mathematical functions that relate the variables and responses of Tables 6 and 7, using response surface methodology.
Taking into consideration the conclusions drawn in the previous section, referring to the collective influence of the design variables on te responses, three different central composite designs have been developed, as follows: Each central composite design has been developed through 2 k points from the factorial design with k factors; 2k axial points face centered, where one variable takes the upper and lower limits and the others have mean values; and finally one central point. Thus, a total number of 65 experiments (15,25, and 25, respectively) have been completed. Also, prior conducting the experiments, the design variables must be normalized to values (−1), (0), and (+1), which stand for the lower bound, mean value, and upper bound of each variable, respectively (Equation (9)).
where k UPi and k LOWi are the upper and lower limits defined in Table 5. Firstly, the study has focused on the relationship between f FEA1 and the variables that affect it. Using the results obtained from the central composite design, an initial second-order model with all of the design variables and interactions have been developed, namely Equation (10), as follows: where f RSM1 is the estimated response corresponding to the first natural frequency, f FEA1 .
Then, the significance and the predictive capability of the regression model as well as the significance of the individual regression coefficients have been examined by means of the coefficients of determination, R 2 and t-tests (Tables 8 and 9). In Table 8, the R 2 coefficients for the initial model show that the regression function explains the observed responses in the central composite design experiment quite well. Also, pred R 2 suggests that the model will fit new responses remarkably.
In order to test the significance of the different terms of Equation (10), the t-statistics [17] for coefficients b j of the initial model have been calculated (Table 9). Using a 95% confidence level (α = 0.05), these terms must be larger than the value of the t-distribution t 0.025,5 = 2.571, and it is shown that the corresponding t-statistics for b 13 , b 23 , and b 33 are smaller. Thus, these three terms are non-significant in the regression model and can be removed. As it is convenient to eliminate one term in each step, b 13 has been picked out first, as its t-statistic was the smallest one.
Then, the same procedure, explained in previous paragraphs, has been repeated for this new model, without the term b 13 . The results in Table 8 (model 2) show that both adj R 2 and pred R 2 have increased slightly. Therefore, as expected, removing the non-significant terms in the regression model has led to a more adequate model. Nevertheless, in the regression equation still there are non-significant terms (Table 9), as some t-statistics are smaller than t 0.025,6 = 2.447. So, coefficient b 23 has been removed, and a new model (model 3) has been made. In this case, adj R 2 and pred R 2 have reduced slightly. Although the differences are totally negligible, considering that this regression model would lead to poorer results than the previous one, the iteration process has been stopped and the preceding regression model has been selected. In Table 10, the results of analysis of variance (ANOVA) of the final model are summarized.  Table 10, it is shown that the model is highly significant (p < 0.001), and confirms that it can be used to simulate the response adequately.
So, the final regression equation for the first natural frequency is as follows: f RSM1 = 34.9265 + 1.1893 · X Y3 + 1.3810 · X Y22 + 0.1806 · X Z4 + 0.1531 · X Y3 · X Y22 + + 0.0177 · X Y22 · X Z4 − 0.4593 · (X Y3 ) 2 − 0.5171 · (X Y22 ) 2 − 0.0559 · (X Z4 ) 2 (11) In this case, only one b j element has been removed and the optimum model is very similar to the initial one. As it will be shown later, in some regression equations, more b j elements will be eliminated, mainly those referred to in second-order terms (see, for example, Equation (17)).
A similar procedure has been followed for f FEA5 . In this case, the regression equation is as follows: f RSM5 = 87.4450 + 0.7551 · X Y9 + 0.5414 · X Y22 + 2.3563 · X Z13 + 0.7345 · X Z63 + + 0.0618 · X Y9 · X Y22 + 0.1096 · X Y9 · X Z63 + 0.0785 · X Y22 · X Z63 + 0.0725 The coefficients of determination, R 2 = 0.9989, adj R 2 = 0.9977, and pred R 2 = 0.9984 (Table 11), have led again to a reliable model. Table 11. Summary of coefficients of determination. RSM-response surface methodology. Finally, the rest of the responses, f FEA3 , f FEA4 , and f FEA6 , along with MAC 3 and MAC 4 , have been analyzed altogether, because the variables that affect them were the same. The final regression equations are shown in Equations (13)- (17), and the coefficients of determination in Table 11.
MAC RSM3 = 63.4984 − 9.4610 · X X8 + 8.9519 · X X210 + 2.1027 · X Y3 − 3.8195 · X Z4 + + 10.4715 · X X8 · X X210 + 1.1365 · X X8 · X Y3 − 2.6332 · X X8 · X Z4 + 3.0198 · (X Z4 ) 2 (16) MAC RSM4 = 84.2937 + 5.9049 · X X8 − 5.2622 · X X210 − 1.0676 · X Z4 + 12.2980 · X X8 · X X210 − − 9.2492 · (X X8 ) 2 (17) In Table 11, the coefficients of determination for f RSM1 , f RSM5 , and f RSM6 are very close to 1.0, while the coefficients for f RSM3 and f RSM4 are slightly lower, although greater than 0.974, and all of them are similar or better than those attained by the authors of [22][23][24]. On the other hand, the coefficients of determination for MAC RSM3 and MAC RSM4 are lower, in some cases under 0.9. In this case, it is not possible to compare them to others, because, to the best of our knowledge, in the literature, there are no results using RSM to simulate MAC responses. Nevertheless, those values are also superior to the coefficients obtained by the authors of [22][23][24] for other responses. So, in conclusion, the approximate functions in Equations (11)- (17) were judged as good enough to accurately relate the design variables and responses, and are adequate to use in the subsequent phase of the improvement procedure.

Determination of Updated Values of Design Variables
Once the explicit relationships between the design variables and model responses have been determined, the next step is to identify the most adequate stiffness values for the connection elements of the FE model, so that the new model simulates accurately the experimental dynamic behavior.
However, prior to performing this step, it is interesting to have a look at Table 12, which shows the responses obtained and the combinations of variables in the CC design 2, and wherein the frequencies inside the range (f expi − 1 Hz) < f FEAi < (f expi + 1 Hz), and where MAC values higher than the initial ones have been bolded. These facts suggest that it is not possible to develop a FE model that fits those three frequencies and that provides acceptable MAC values with a unique value of design variable k X8 . In fact, Wu et al. [36] have also addressed a similar behavior in other machine tool with roller type linear guideways. Therefore, it will be necessary to identify one k X8 value to match, in combination with k X210 , k Y3 , and k Z4 ; natural frequencies f FEA3 , f FEA6 ; and necessarily MAC 3 , as well as other k X8 value to match f FEA4 and MAC 4 , taking into account that design variable k X8 does not affect the rest of the responses ( Table 6).
For that purpose, the desirability function has been used, as explained in Section 3.3. Two types of desirability functions have been defined (Figure 9), one for natural frequencies and the other one for MAC values.
where MAC0 has been selected taking into consideration Tables 4 and 12.
Finally, the global desirability function D, Equation (22), is composed by combining one global desirability function for the frequencies, Df, and another function for the MAC values, DM, and applying weighting coefficients wf and wM to each of them. Table 13 shows the updated values of the normalized design variables Xi, and the corresponding natural values ki. As mentioned before, two stiffness values kX8 have been estimated, namely: (1) is valid for the frequency ranges from 0 Hz to 72 Hz, and from 100 Hz to the upper limit of the range of interest, and (2) is adequate for the remaining range, which includes the 4th natural frequency.  These functions can be expressed in mathematical form as follows: where MAC 0 has been selected taking into consideration Tables 4 and 12.
Finally, the global desirability function D, Equation (22), is composed by combining one global desirability function for the frequencies, D f , and another function for the MAC values, D M , and applying weighting coefficients w f and w M to each of them. Table 13 shows the updated values of the normalized design variables X i , and the corresponding natural values k i . As mentioned before, two stiffness values k X8 have been estimated, namely: (1) is valid for the frequency ranges from 0 Hz to 72 Hz, and from 100 Hz to the upper limit of the range of interest, and (2) is adequate for the remaining range, which includes the 4th natural frequency. These values have been driven into the FE model and the posterior FE analysis has led to the frequencies and MAC values, indicated in Table 14. For the sake of comparison, the simulated responses, f RSMi and MAC RSMi , obtained in Equations (11)-(17), when the updated values of the design variables are substituted, are also shown. From Table 14, it can be seen that the quadratic regression equations have provided values of the simulated frequencies that almost coincide with the values obtained after the completion of the FE analysis. In fact, the maximum distance is 0.4 Hz in the 6th frequency, which is really insignificant. The difference in the simulated MAC values is greater, which is in accordance with Table 11, where it was suggested that the predictive capability of the MAC regression equations was inferior. Therefore, both the regression meta-models and also the statistic indicators, R 2 and t-statistic, have performed adequately.
Finally, once the identified values of design variables have been incorporated into the FE model, the resultant dynamic responses have shown a closer match to the experimental results, proving the adequacy of the conducted procedure. Thus, in Figure 10, two synthesized FRFs obtained from the updated FE model are compared to the corresponding experimental FRFs in reference point 5 (Figure 2), and the agreement is quite reasonable.
From Table 14, it can be seen that the quadratic regression equations have provided values of the simulated frequencies that almost coincide with the values obtained after the completion of the FE analysis. In fact, the maximum distance is 0.4 Hz in the 6th frequency, which is really insignificant. The difference in the simulated MAC values is greater, which is in accordance with Table 11, where it was suggested that the predictive capability of the MAC regression equations was inferior. Therefore, both the regression meta-models and also the statistic indicators, R 2 and t-statistic, have performed adequately.
Finally, once the identified values of design variables have been incorporated into the FE model, the resultant dynamic responses have shown a closer match to the experimental results, proving the adequacy of the conducted procedure. Thus, in Figure 10, two synthesized FRFs obtained from the updated FE model are compared to the corresponding experimental FRFs in reference point 5 (Figure 2), and the agreement is quite reasonable.

Conclusions
This paper presents a consistent methodology to improve the FE models of complex mechanical systems using, in an integrated way, different numerical and experimental techniques. The procedure is applied in a machining center with numerous uncertainties in the internal connections and supporting conditions. In this methodology, the complete design space is encompassed, so that the selection of the initial values of the design variables, which is one of the major drawbacks of the sensitivity-based methods, due to the variable sensitivity values along the design space, is avoided.

Conclusions
This paper presents a consistent methodology to improve the FE models of complex mechanical systems using, in an integrated way, different numerical and experimental techniques. The procedure is applied in a machining center with numerous uncertainties in the internal connections and supporting conditions. In this methodology, the complete design space is encompassed, so that the selection of the initial values of the design variables, which is one of the major drawbacks of the sensitivity-based methods, due to the variable sensitivity values along the design space, is avoided.
Firstly, it is demonstrated that the two-level fractional factorial design is an effective tool to perform parameter screening, as the most significant parameters and two-factor interactions are detected. For this purpose, instead of using one cumbersome resolution V design, more simple fractional designs with fewer parameters are gradually completed and examined. This procedure allows for removing high-order interactions and to circumvent the presence of constraints between the variables, and leads to a better comprehension of the influence of the design variables on the behavior of the mechanical system. Also, in this step, it is shown that the design variables perform a kind of collective work, as groups of them affect groups of responses. In addition, some of them can be satisfactorily removed, in contrast to other findings reported in literature, where it is assumed that the complete set of the selected design variables is significant. This is a key feature of the proposed methodology, as it allows for diminishing the complexity of the subsequent regression equations, due to a substantial drop in the number of terms.
In the second step, it is demonstrated that the relationships between the stiffness parameters and the modal responses of the machine tool can be accurately expressed by second-order functions. A combined procedure using coefficients of determination and the t-statistic is applied to remove the non-significant terms, and thus the accuracy of regression equations is increased. At this point, the assessment of the predictive capability of the regression meta-models plays an important role. The regression equations for the MAC values are also used to ensure the correspondence between the numerical and experimental responses, because if only the frequency values are matched, it could lead to unacceptable MAC values. So far, this issue has been overlooked in the literature.
Also, the use of central composite designs allows for developing quadratic regression equations at an acceptable cost. In addition, these designs have led to detecting that the stiffness of one connection is dependent on the relative movement between the modules of the machining center.
It is proved that the quadratic regression equations are adequate to accurately identify improved values of the design variables, because when these values are included in the FE model, minimal differences between the FE and experimental responses are found. Also, because of the substitution of the full FE model by polynomial functions, the identification step, usually a costly iterative procedure, is accelerated. The use of desirability functions and weighting factors facilitates the progress of this step.
Finally, the presented methodology can be generalized to any machine tool and any design variable (damping in connections, Young's modulus, mass density, etc.) and will allow for obtaining an updated finite element model, which would serve as a starting point to optimize the machine design and eliminate stability problems under operating conditions. Potential future research directions include the analysis and implementation of other designs (orthogonal, Latin Hypercube, and D-optimal) for fitting second-order models with constraints in the design space. Another possible direction is the use of techniques of model reduction or transformation to modal space to diminish the time-consuming DoE runs. Funding: This work has been supported by the University of the Basque Country UPV/EHU under programs PPGA17/04 and US17/16, and the Spanish Ministry of Economy, Industry, and Competitiveness under program DPI2016-74845-R.

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