Learning-Based Prediction of Pose-Dependent Dynamics

: The constantly increasing demand for both, higher production output and more complex product geometries, which can only be achieved using ﬁve-axis milling processes, requires elaborated analysis approaches to optimize the regarded process. This is especially necessary when the used tool is susceptible to vibrations, which can deteriorate the quality of the machined workpiece surface. The prediction of tool vibrations based on the used NC path and process conﬁguration can be achieved by, e.g., applying geometric physically-based process simulation systems prior to the machining process. However, recent research showed that the dynamic behavior of the system, consisting of the machine tool, the spindle, and the milling tool, can change signiﬁcantly when using different inclination angles to realize certain machined workpiece shapes. Intermediate dynamic properties have to be interpolated based on measurements due to the impracticality of measuring the frequency response functions for each position and inclination angle that are used along the NC path. This paper presents a learning-based approach to predict the frequency response function for a given pose of the tool center point.


Introduction
Chatter vibrations are common challenges in milling processes, leading to an insufficient workpiece quality and reduced lifetime of the machine tool and cutting tools [1,2], especially if long and slender milling tools are necessary to machine the desired workpiece geometry [3][4][5][6]. Several approaches can be used to optimize the dynamic behavior of milling operations. The process stability can be evaluated analytically based on measured or simulated frequency response functions (FRFs) [7][8][9]. Furthermore, simulation approaches offer the possibility to reduce run-in periods of identifying suitable process parameter values and, thus, accelerate and simplify the process design and optimization, even for high process run-times [2,10,11]. For an optimization of milling processes with varying engagement conditions and complex desired workpiece shapes, geometric physically-based process simulations can be used [10]. In this context, the dynamic behavior of the compliant system, consisting of the combination of the machine tool, spindle, and cutting tool, can be modeled by a set of uncoupled, damped harmonic oscillators to represent the FRF of the system measured at the tool center point (TCP) [12]. For each oscillator of the set, the modal mass m m , the natural frequency f m and the damping constant γ m have to be identified. Subsequently, the process forces, exciting the compliance system and estimated by, e.g., applying an empirical force model, which has to be calibrated based on the used combination of tool geometry and workpiece material, can be used within the dynamic model. Thus, the deflections of the milling tool can be calculated and the process stability can be assessed. Additionally, stability lobe diagrams (SLD) can be calculated for different combinations of the process parameter values using suitable stability criteria. However, especially when machining free-formed surfaces of large workpieces, the pose-dependent load of the spindle bearings and axis drives influences the modal properties of the system significantly [11,[13][14][15], resulting in different frequency response behaviors for each pose defined by the NC path.
Different methodologies have been investigated to model this influence. Kono et al. [16] investigated the influence of the rotation angle of a five-axis machine tool with swiveling table on the dynamic compliance and oscillations during a milling process by modal analysis. The maximum difference between the compliance amplitudes was about 40 % when varying the B-axis angle. Multiple analytical approaches can be found in literature. In this regard, Budak et al. [17] presented analytical models for the cutting geometry, process forces and process stability for 5-axis milling operations. Shamoto and Akazawa [18] analytically calculated SLDs based on FRF measurements for a ball-end milling process with respect to different inclination angles. A reasonable agreement between measured and calculated stabilities could be achieved. When considering the modeling of pose-dependent dynamics, Du et al. [19,20] modeled the dynamic behavior of a bi-rotary milling head using multi-rigid-body dynamic models considering the varying stiffness of the flexible joints affected by gravity and cutting forces. Furthermore, Chao and Altintas [21] optimized the tool path for machining a free-formed surface by analytically analyzing the engagement conditions using FRFs which were measured in different poses. Regarding oscillations of thin-walled workpieces, Siebrecht et al. [22] used a barycentric interpolation approach for estimating the dynamic compliance of the workpiece at different positions on the surface.
In recent years, the use of machine learning (ML) methods in research activities regarding production engineering has increased drastically [23][24][25]. They offer the capability to deliver predictions of process characteristics for previously unseen process features and conditions with a reasonable accuracy. Thus, ML-based models were successfully trained in order to predict process forces [26][27][28][29], the surface roughness of the machined workpiece [30][31][32] and tool vibrations [33,34]. When considering SLDs, Friedrich et al. [35] realized an ML-based framework that was able to iteratively adapt an SLD based on continuous measurements during process conduction. In addition, Denkena et al. [36] learned SLDs based on measured data and investigated the suitability of different methods. There are few publications dealing with the prediction of pose-dependent dynamics based on ML methods. In this context, a transfer learning approach was presented by Chen et al. [37], whereby a multilayer perceptron (MLP), which was trained using a high amount of data originating from impact hammer tests, was adapted to be valid for different cutting tools using only few additional data. However, FRFs with only one single dominant mode were considered and no SLDs were derived for the regarded configurations.
In this paper, an investigation of the prediction of FRFs and the subsequent fitting of corresponding oscillator-based compliance models is presented in order to enable the interpolation between poses for which FRFs were acquired by impact hammer tests. Using an evolutionary-based algorithm, an automated parameterization of compliance models using an initial set of parameter values, estimated by a single manual fitting procedure, could be derived to reduce the amount of manual fitting procedures for given FRFs. Furthermore, SLDs were generated using the resulting compliance models and a geometric physically-based simulation approach. These SLDs are compared to the experimentally determined stability limits using acoustic emission signals for different inclination angles, visualizing the usefulness of the approach for an optimization of five-axis milling operations. The SLDs are also compared to SLDs that were acquired using compliance models predicted by ML models directly, skipping the investigation of FRFs at the interpolation poses. In addition, the advantages and disadvantages of both approaches are discussed in detail. Different methods of ML were used for the learning objectives to evaluate their suitability. The paper is structured, as follows. In Section 2, the technological investigations, which were conducted to acquire FRFs for different poses using two machine tools and the acoustic emission signals for different inclination angles, are described in detail. The geometric physically-based simulation system, which was applied to estimate SLDs, is introduced in Section 3. The identification of the modal parameter values of the oscillator-based compliance models to represent FRFs while using an evolutionary algorithm is explained in Section 4. Section 5 presents the learning objectives and the used base-line methods, which were applied to perform the interpolation of FRFs and compliance models. The results of the learning tasks for two different five-axis machine tools are presented, evaluated and discussed in Section 6. In addition, Section 6 also comprises the results of comparing SLDs for predicted FRFs as well as predicted compliance models to stabilities, retrieved by an evaluation of acoustic signals of corresponding milling processes. The paper concludes with a summary of the conducted investigations (see Section 7).

Technological Investigation
For the acquisition of a sufficient training set and, therewith, validation of the presented method, frequency response measurements were performed through impact hammer tests on two machine tools in different poses. On a five-axis machining center with swivel head kinematic Heller FT 4000, denoted as M 1 in the following, the positions of the axes influencing the pose of the tool, X, Y, and C, were varied. Furthermore, on a five-axis machining center with fork head kinematic DMG HSC 75 linear, denoted as M 2 , the Y-, Z-, and B-axis positions were varied. For both of the machining centers, the axis positions were selected on the basis of an extended centrally composed experimental design with star points, in order to be able to investigate higher order influences as well as possible cross-correlations (cf. Figure 1). A total number of P 1 = 46 and P 2 = 49 poses were investigated using M 1 and M 2 , respectively. All of the axis positions were defined in the machine coordinate system allowing the axes to be moved individually over the whole working area of the machine tools (cf. Table 1). Range of investigated axis positions X −400 mm to +400 mm Y 0 mm to +600 mm Y +100 mm to +900 mm Z −400 mm to 0 mm The FRF measurements were conducted using a Fraisa X7400 spherical end mill with a diameter of d = 10 mm and a free length of l l = 30.4 mm. Hereto, an impact hammer Kistler 8206 and an acceleration sensor 352C23 by PCB Piezotronics, mounted at the tool tip, were used. Fast Fourier transformations were performed up to f max, 1 = 6400 Hz and f max, 2 = 3200 Hz for M 1 and M 2 , respectively. The finite length of the impulse response function measurements led to a frequency resolution of ∆ f 1 = 0.25 Hz and ∆ f 2 = 0.5 Hz for M 1 and M 2 , respectively.
Slot milling experiments were conducted using the aforementioned milling tool, the machine tool M 2 , and an increasing depth of cut a p , using a starting value of a p,start and a depth of cut of a p,end at the end of each slot, in order to validate the calculated SLDs. A high speed steel ASP 2012 hardened to approx. 58 HRC was machined. Different inclination angles and spindle speeds were investigated. There are different publications in literature, which investigated the detection of chatter vibrations applying wavelet analysis using various signals, e.g., cutting forces [38] or tool accelerations [39][40][41]. In this contribution, discrete acoustic emission signals were recorded and analyzed utilizing the continuous wavelet transform [42] for each experiment. In this context, ξ (i) is the i-th sample of the original acoustic emission signal, N is the total number of samples of the signal, Ψ * is the complex conjugate of a mother wavelet Ψ, δ t is the time difference between two samples of the signal, and a and b are the scaling and translation variables, respectively. Each scaled and translated mother wavelet corresponds to an investigated frequency at a time instant. By calculating the convolution between the original signal and the scaled and translated mother wavelet, the correlation between the signal and a frequency can be estimated for each time instant. For the mother wavelet, the complex Morlet wavelet [43] Ψ(η) = π − 1 4 e iω 0 η e − η 2 2 (2) was used. The stability limit a p,crit of each experiment was estimated as a p,crit = a p,start + i crit · a p,end − a p,start whereby W(a) is the wavelet transform of the signal, τ is a user-defined threshold value and S is the set of investigated scales, which correspond to a set of natural frequencies of the FRF of the regarded pose. For the Morlet wavelet, the relationship between a scale a and the Fourier period λ can be described as whereby ω 0 is the non-dimensional central frequency of the mother wavelet. Because natural frequencies measured for a non-rotating spindle may differ from natural frequencies of a rotating system, the function T(W (i) (a), N T ) calculates a weighted average of the wavelet intensities, incorporating N T neighboring frequencies of the natural frequency, which corresponds to the scaling a.
For the weighting, the Blackman window function [44] w(j) = α 0 − α 1 cos 2πj with α 0 = 0.42, α 1 = 0.5, α 2 = 0.08 (7) was used, whereby W (i) (a) was weighted with w(N T /2). Equation (4) represents different values for the translation variable b with the index i, assuming that there is a set W (i) (a) ∀i ∈ {1, . . . , N} , whereby each W (i) (a) corresponds to a sample ξ (i) of the original signal and b is estimated according to the temporal location of ξ (i) within the time series. Thus, the stability limit is estimatednby identifying a critical index i crit , which is defined as the index, where each of the weighted averages of wavelet intensities of a frequency and its neighboring frequencies of a set of investigated natural frequencies exceeds W(a) + τ · W s (a), whereby are the mean and the standard deviation of a wavelet transform W(a), respectively. Figure    Acoustic emission signals often are highly influenced by noise. As a result, a high amount of frequencies correlate with the corresponding scaled and translated wavelet transforms. The natural frequencies 1415 Hz, 1500 Hz and 1580 Hz of the three dominant peaks of the FRF in X-and Y-direction were considered when analyzing the wavelet transform, since they were expected to have the most influence on the stability behavior of the process. Using the aforementioned approach, a stability limit of a p,crit = 0.5 mm was identified, since all og the wavelet transforms, which corresponded to the three natural frequencies, exceeded the defined threshold at the corresponding point in time. The stability limits, which are calculated using the presented method, have to be interpreted with caution. The resulting limits are highly dependent on the signal-to-noise ratio of the acoustic emission signal and the choice of τ. As a result, the transferability to processes that differ from the considered configuration is limited.

Geometric Physically-Based Simulation of Milling Processes
In general, geometric physically-based process simulations combine geometric representations of the milling tool and workpiece with physically-based models in a time-discrete manner in order to offer predictions for various process characteristics [2]. These characteristics usually comprise macroscopic effects such as process forces, process dynamics and the resulting surface errors of the machined workpiece surface, which can be modeled by incorporating and combining different models [10,45]. Apart from milling processes, grinding [46] and honing [47] processes or process chains [48] can also be analyzed.
In this contribution, a simulation system was used, which is described in detail by Wiederkehr and Siebrecht [2]. The constructive solid geometry (CSG)-technique [49] was used to analyze the uncut chip shape for each simulation time step by intersecting the models of tool and workpiece (cf. Figure 3) [50]. The resulting process forces were calculated using an empirical force model [51,52], which can be described as for the cutting, normal and tangential direction, respectively, whereby h 0 = 1 mm, h is the uncut chip thickness, b is the width of the cutting slices, which represent the spacial discretization of the cutting edges of the geometric model of the cutting tool and k c , k n , k t , c c , c n , and c t are the model parameters.
Machining experiments were conducted using the process parameter values defined in Table 2 in order to parametrize the model parameter values. The force model parameter values were identified using a gradient-based optimization algorithm [53]. Because the uncut chip thickness remains constant for different inclination angles and its nonlinear influence on the process forces is insignificant for the subject under investigation a linear cutting force model could be assumed. Thus, the exponent in Equation (10) was set to 1, i.e., c i = 0.0. The identified parameter values of the force model for different inclination angles are summarized in Table 3. Table 3. Parameter values of cutting force model. A compliance model based on uncoupled, damped harmonic oscillators for representing FRFs was used for simulating the resulting deflections of the tool and evaluating the process stabilities [12,54]. As discussed by Surmann and Enk [54], a set of Poincaré points were used to calculate the diameter of the Poincaré section, which was used as stability criterion. Each Poincaré point was defined two-dimensionally as the magnitudes of the deflections in Xand Y-direction at the beginning of each tooth engagement. The diameter of the Poincaré section was then calculated as the diagonal length of the smallest rectangle, which could be defined around the set of considered Poincaré points. The smaller the diameter, the higher the assumed process stability. For each combination of depth of cut and spindle speed in the regarded range of the parameter values, simulations were conducted and the resulting stabilities were evaluated in order to generate SLDs for a given inclination angle.

Fitting of Parameter Values of Compliance Models
In order to represent FRFs, either measured or predicted compliance models were used in Xand Y-direction of the tool coordinate system according to the method proposed by Surmann et al. [12]. These consisted of a set of uncoupled damped harmonic oscillators separately for the Xand Y-direction. Each oscillator was parameterized by identifying values for the modal mass m m , the natural frequency f m , and the damping constant γ m . The parameterized models were used to calculate the complex response function (11) according to the amplitude and phase resulting from each oscillator q for each investigated angular frequency ω. The results were then compared to the measured FRFs to calculate the optimization loss. Therefore, the resulting problem is an optimization problem of parameter values and can be solved by a genetic algorithm [55]. Genetic algorithms are inspired by the concept of natural evolution and use selection, mutation, and crossover methods to evolve a population of solution candidates to retrieve an optimized result [56,57]. The candidates are called individuals and consist of a set of parameter values which is a possible solution for the problem, where the quality is defined by a fitness function [56]. For the oscillator fitting, one individual is a set of 3 · Q parameters, where Q is the number of oscillators. To initialize the individuals of the starting population of the automated fitting procedure using the genetic algorithm, oscillators for the Xand Y-direction were generated manually for one FRF of each machine. Subsequently, the fitness values were calculated according to the inverse of which represents the sum of the squared error between the normalized values of the calculated amplitudeÃ p and phase valuesφ p and the normalized measured dataÃ t andφ t for the investigated frequency range. The normalization was performed according to the minimum and maximum values of the amplitudes and phases to transform the values, which were used by the fitness function, to the range 0-1. Subsequently, the individuals of one population were sorted according to the corresponding fitness values in ascending order. Operations of selection, mutation, and crossover were used to estimate the next generation of individuals. The mutation function added a random value between −2 Hz and 2 Hz to f m of the considered oscillator and used a deviation of up to 2 % to modify m m and γ m . Applying the crossover method, the oscillators of two selected individuals, from 30 % of the actual generation that achieved the best fitness values, were mixed. This procedure was repeated until a total amount of 1000 iterations were conducted or if the best fitness value did not improve for 20 iterations. An analysis of the measured FRFs resulted in the assumption that there were frequencies above the considered frequency measurement range, which were relevant to achieving a reasonable fitting result at the upper frequency limit. Therefore, an additional oscillator with fixed parameter values was added, whose natural frequency was above the measurement range. Figure 4 shows an exemplary comparison of the measured and calculated FRFs, where the calculation was based on an evolutionary parameter optimization of oscillator-based compliance models.

Measurement range
Measurement range (a) (b) Figure 4. Comparison of measured and calculated FRFs, where the calculation is either based on compliance models whose natural frequencies (a) entirely lie inside the considered frequency range or (b) use exactly one oscillator that had a fixed natural frequency outside the measured upper frequency limit.
The calculated FRFs in Figure 4a are based on compliance models whose natural frequencies lie inside the considered frequency measurement range. Unfortunately, using this strategy, the fitting results were not suitable to accurately represent the measured frequency response behavior. There clearly was a phase angle shift at the end of the frequency range, indicating the necessity of another oscillator whose natural frequency is outside of the frequency range. Because no reference measurements were present to consider this oscillator within the parameter optimization, an oscillator with the fixed parameter values f m = 3807.51 Hz, γ m = 2236.61 1 s and m m = 0.04 kg for X-direction and f m = 3582. 16 Hz, γ m = 1527.42 1 s and m m = 0.09 kg for Y-direction was included to the oscillator set of each compliance model. Using this approach, a suitable fit could be achieved for each FRF of each considered pose. Nevertheless, there was an uncertainty regarding the accuracy of the representation of the frequency response behavior above the frequency measurement range, which was neglected in the following investigations.

ML Methods for Predicting Pose-Dependent Dynamics of Milling Processes
Two learning objectives were considered in this paper. For both objectives, let X be a set of J × N features, sampled from an unknown distribution D, and Y be a set of K × N targets, labeled by some target function. The first objective comprised the prediction of FRFs. For this, the measured FRFs of all P considered measuring poses were discretized into data points by the frequency resolution ∆ f , so that M is the number of investigated frequencies for each pose. Each of the N = P · M data points contained a number of K targets, comprising the compliance amplitude and phase shift for the Xand Y-direction of the machine coordinate system. Let J be the number of features, consisting of the frequency and positions of the three axes that define the pose. For the second learning objective, which represented the prediction of modal parameter values for given poses, let J consist of the three pose-dependent features and N = P. For the targets, let K = 3 · (Q x + Q y ), whereby Q x and Q y are the number of oscillators in Xand Y-direction, respectively. Using this approach, the learning task tried to also represent the relationship between different interdependent oscillators of each compliance model, even across the two different oscillation directions. For both learning objectives, the goal was to find a learner h : X → Y with respect to the distribution D. Because D is unknown to the learner, the true prediction error cannot be calculated. To mitigate this problem, the empirical risk L D (h) was employed. In practice, it is desirable to find a learner, such that is minimized. This paradigm is often denoted as empirical risk minimization. Several learning-based methods can be used for approaching regression objectives. The simplest models can be trained while using linear regression [58,59], whereby the relationship between the features and the targets are assumed to be The parameters β are the coefficient of the model and the error term represents all of the deviations between forecasts and Y which cannot be represented by the model. One of the most popular approaches for fitting a linear model to a set of training data is the least-squares method. In this context, the approach tries to identify the model coefficients β, which minimize the sum of squared residuals Setting the first derivative of Equation (15) to zero, the unique solution is given bŷ if X T X is invertible. The simplicity of linear regression, which offers highly interpretable model instances, entails several issues, e.g., either delivering a poor prediction accuracy or the trained model suffers from severe overfitting. This bias-variance dilemma [60] can be addressed through several regularization techniques. In this contribution, the elastic net [61] was used as a regularization approach, which combines the 1 and 2 norms of the model coefficients, such that where λ 1 and λ 2 are regularization parameters. Using this approach, the elastic net tries to overcome issues, which arise when using either ridge [62] or least absolute shrinkage and selection operator (LASSO) [63] regression. While penalizing the size of β through the 2 norm in the context of ridge regression leads to lower variance, coefficients of highly correlated variables tend to be shrunk together.
Using LASSO, a sparse solution is encouraged by incorporating the 1 norm, which aims to make the model more interpretable. However, although no shrinkage of coefficients is done, LASSO simply picks one of the group of coefficients of highly correlated variables. No systematics of the selection can be derived, depending on the context. Elastic net attempts to stabilize this selection procedure by shrinking and selecting sparsely simultaneously, realizing a compromise between ridge and LASSO. Another popular approach to reduce overfitting effects is the use of ensemble techniques [64,65], whereby a set of weak learners h b (x) with high bias are utilized to generate a strong combined learner. Usually the weak learners are tree-based models mostly trained using the classification and regression trees (CART) algorithm [66], but the choice of the used model is arbitrary. Several strategies emerged for the combination of the weak learners. In bootstrap aggregation (bagging) [67], each weak learner is trained using a different subset of the training set. This subset is sampled with replacement, which means that an observation can be in multiple subsets. For regression objectives, the final prediction of the combined learner can be summarized aŝ In this paper, random forests (RF) [68] were used, which can be categorized to bagging techniques. While training a RF, opposing to the algorithms of classical decision tree learning, the split is performed using a randomly sampled set of features, also known as feature subsampling. Using boosting [69], a sequential approach is performed for the combination of the weak learners. At each iteration, the model tries to learn from mistakes of the previous iteration, whereby the step length ρ b is usually estimated by line search. As a popular example of boosting algorithms, XGBoost [70] was additionally used in this contribution. XGBoost is a specific parallelized implementation of the gradient boosting (GB) technique [71], whereby each weak learner is fitted to the residuals of the previous learner by formulating the objective as a gradient descent optimization procedure. The regularized objective , penalizes both the number of leaves T and the size of the weights β of the leaves of each tree h b , whereby l is a loss function, representing the difference betweenŷ i and y i . This addition aims to reduce overfitting in a similar manner than using the ridge regularization in the context of linear regression. Furthermore, feature subsampling, originating from RF, is used when performing the tree split to further reduce overfitting and speed up the training procedure. An approximation algorithm is conducted to improve the runtime efficiency when utilizing tree boosting using the XGBoost algorithm, instead of performing an exact greedy search, which iterates over all possible splits over all features. This algorithm tries to find the best split that identifies a set of candidate splits according to the proposed distributed weighted quantile sketch algorithm.

Results
The following section presents the results for the two considered learning tasks, i.e., the prediction of FRFs and the prediction of parameter values of compliance models, in order to reduce the measurement effort and to enable the possibility to retrieve information about the pose-dependent dynamic behavior of milling operations for poses that were not investigated technologically. For both learning tasks, different strategies for defining the training and test sets were investigated. As described in Section 2, FRFs were measured for a total amount of 46 and 49 poses while using M 1 and M 2 , respectively. In addition to randomly choosing four poses for testing, another strategy for choosing poses for the test set was to use all poses that use a certain inclination angle, whereby three different test angles were considered in total for each machine and each investigated base-line method. This resulted in a total amount of five test poses for each test angle of each machine tool, except for a test angle of B = 45 • for M 2 , where four test poses were examined. In order to measure the deviation between measured and predicted data, the root-mean-square error (RMSE) for each target y k was used. The results that were achieved by applying the ML base-line methods were also compared to a naïve linear interpolation approach for both learning objectives. In this contest, the pose-dependent features were organized in a k-d tree [72]. For each test pose, a previously specified number of nearest neighboring poses were identified in the k-d tree. Subsequently, each calculated target value for the test pose resulted from a weighted average of the target values from the neighboring training poses. The weights were chosen according to the distance between the test pose and the respective training poses. A parameter that significantly influenced the interpolation accuracy was the number of incorporated neighboring poses. Each possible number of neighbors from a range of 2-20 was examined and the best performing configuration was used to generate the reported results.

Learning of FRFs
In the context of this learning objective, J = 4, since the three pose properties and the frequency were considered to be features. Furthermore, K = 4, comprising the values for the amplitude and phase for the X-and Y-direction. Figure 5 shows a comparison between FRFs that were measured using the machine tool M 1 and predicted FRFs using RF for four different test poses that were not included in the training set. X,Y,C between 2200-2500 Hz. This could be explained by the measured behavior of the phase being quite unique among all measured phases in this frequency range, resulting in the model failing to achieve an accurate prediction. The amplitudes in X-direction were predicted with a nearly non-visible deviation from the measured curves for all test poses. There are two peaks visible in measured FRFs in Y-direction in a frequency range of 1200-1700 Hz, which are very distinct for P (2) X,Y,C and P X,Y,C , the second peak can hardly be detected and for P (4) X,Y,C , the second peak is not present. This behavior could be represented by the model to a certain extent. The fusion of the two peaks for P (4) X,Y,C could successfully be predicted. In addition, the model also predicted a visible second peak for the remaining test poses, but the differences in the distinction between the two peaks across the poses could not be achieved.
This effect was very minor in the investigated data and has to be analyzed in more detail in future research activities. In order to represent such behavior, more observations in which the fusion of peaks is present and that can be considered for the training procedure would be necessary. Figure 6 shows a comparison between measured FRFs using the machine tool M 2 and predicted FRFs using XGBoost, which was the best performing method in this case, for two different test poses. There was a high accordance between the measured and predicted phase angles in Xand Y-direction. The difference between measured and predicted amplitudes in Y-direction was also very low. Although the width and hight of the amplitudes in X-direction could successfully be predicted, there was a deviation of the shape of the curves in the range between 1300-1700 Hz. This deviation also corresponds to a fusion of peaks, since there are two peaks visible for P (1) Z,Y,B , which could not be accurately predicted, and only one peak for P (1) Z,Y,B , whereby the difference between measurement and prediction is low. This emphasizes the already concluded statement about more observations being necessary to predict these effects. Table 4 shows the results of learning FRFs for the two investigated machine tools while using different methods and different inclination angles for testing, which were not used during training.
In general, using the elastic net regression was not successful for both machine tools, since the RMSE is about 10 times higher than using other methods. Surprisingly, the naïve linear interpolation performed significantly better than the elastic net regression, which is also based on a linear model, although these results were still, on average, approximately 25 % worse than the results of the ensemble methods. This emphasizes the possibility to model the relationship between the axis poses and frequency-based machine tool dynamics linearly with reasonable accuracy. The low prediction accuracy of the elastic net could be a result of an unsuccessful regularization leading to overfitting. When comparing the results of RF with the ones achieved by using XGBoost, the bagging technique performed slightly better on average, but, overall, both ensemble methods delivered a reasonably low RMSE. Predicting FRFs for fixed inclination angles delivered slightly better results than randomly chosen angles. This might be counter-intuitive at first sight. An explanation for this could be the randomly chosen test set comprising scenarios, whose relationship between poses, frequencies, amplitudes, and phases were not represented by the training set. Table 4. Comparison between the root-mean-square error (RMSE) of predicted amplitudes and phases in X-and Y-direction of the machine coordinate system for two machine tools using different methods.

Method
Machine Tool

Learning of Oscillator Parameter Values
In contrast to learning FRF, the learning of parameter values of oscillator-based compliance models, whose results are presented in the following, tried to represent the relationship between information about specific poses and compliance model instances, which would result from parameter optimization procedures. A major advantage of this approach is that it enables the possibility to retrieve information about pose-dependent dynamic behaviors through an evaluation of learned models, which, in contrast to conducting a fitting procedure of oscillator parameter values, can be performed in real-time. In the context of this learning task, J = 3, only comprising the pose-related values. No frequency information has to be provided, since this learning task skipped the prediction of frequency response behaviors and aimed to represent the properties of compliance models directly, given a certain pose. Furthermore, as described in Section 5, K = 3 · (Q x + Q y ), representing the parameter values of a set of oscillators. A number of Q x = 5 and Q y = 4 were chosen to model the measured FRFs by the evolutionary-based fitting procedure, which was necessary to generate the training and test sets, in the Xand Y-direction, respectively. Figure 7 shows an exemplary comparison between a measured FRF, an FRF calculated using a compliance model that results from the fitting procedure, which is described in Section 4, and an FRF resulting from predicted parameter values of a compliance model according to the second learning objective for the pose P Z,Y,B = (−50.0, 500.0, 0.0) and the machine tool M 2 .
Measured data Data of fitted oscillators Data of predicted oscillators X-direction Y-direction Generally, a high accordance was observed. Examining the zoomed-in areas of the FRFs, it can be seen that the shape of the measured FRF could only be coarsely represented by the evolutionary-based fitting procedure. Because the fitted oscillator parameter values served as target values for the learning objective, the FRF, which was calculated based on predicted oscillator parameter values, could not reproduce the measured behavior in higher detail than the FRF, which resulted from the fitting procedure. Nevertheless, there were only small deviations between the fitted and predicted data, which was also concluded based on the results of Table 5. Therefore, the learning of oscillator parameter values directly from given poses could be interpreted as successful. Table 5 shows the results of learning oscillator parameter values directly for given poses for different inclination angles for the test set using different methods.
It can be seen that even the elastic net delivered suitably low RMSE scores, especially for predicting the natural frequency of the oscillator models. The naïve linear interpolation provided the lowest performance and achieved on average approximately 31 % worse results than the results of the ensemble methods. For the damping constant and the modal mass, a superior performance of the ensemble models could be observed, indicating a more non-linear relationship between the pose and the damping constant and the modal mass than between the pose and the natural frequency. Similar to the first learning objective, randomly chosen test scenarios were significantly harder to predict than using a fixed inclination angle for the test set. A trade-off has to be considered summarizing the results of evaluating the two regarded learning objectives. Applying the first objective of learning FRFs, the measured behavior could mostly be represented in detail while deviations between measured and predicted data were observed when a fusion of amplitude peaks of the frequency response was present. A higher prediction accuracy should be expected when acquiring a higher number of observations for the training set, which more prominently represent the behavior of fused peaks. If an application of a geometric physically-based simulation system is pursued to, e.g., simulate SLDs or visualize simulated location errors on the workpiece surface based on an evaluation of oscillator-based compliance models, a fitting of the parameter values of the compliance model has to be performed following each prediction of an FRF. Because this has to be conducted for each pose of the used NC path, the run-time of simulation executions would increase significantly, negating one of the major advantages of geometric simulation approaches. Nevertheless, this approach could be very useful for applications in methods that perform stability calculations in the frequency-domain. These are, e.g., the analytical stability calculation method by Altintaş and Budak [7] or the impulse dynamic subspace method proposed by Dombovari et al. [73] used to calculate time-domain-based process simulations directly based on frequency response functions. However, predicting parameter values of oscillators directly without considering the underlying FRF given a pose of the NC path by realizing the second learning objective, no deterioration of simulation run-times is expected. This is due to the capability to evaluate ML models in real-time, so that no fitting procedure is required. By using this approach, a more coarse representation of the frequency response behavior has to be tolerated that results from the automatic evolutionary-based fitting procedure, which was used to generate the training set for the learning task. The fitting procedure has to be performed in both cases, but at various stages of the realization sequence, i.e., following the ML-based prediction of each FRF for the first learning objective or following the measurement of each FRF to generate the training set in the context of the second learning task. Thus, the loss of information is inevitable, relativizing the disadvantages of the second learning objective to a certain extent. However, information of the frequency response behavior for a given pose, where no FRF information was measured, can only be acquired by using the first learning objective. Finally, even a naïve linear interpolation of target values for both learning objectives delivered acceptable results, although it performed worse overall than the ML-based methods. If a time-consuming training procedure of model instances is to be avoided at all costs, this is still a valid approach for generating usable results. Figure 8 shows a comparison between SLDs that were calculated by using a geometric physically-based simulation system (see Section 3), fitted compliance models (see Section 4), resulting from either using measured (see Section 2) or predicted FRFs and predicted compliance models (see Section 5). In addition, the SLDs are compared to experimentally identified stability limits according to the approach that is described in Section 2. In this context, τ ∈ [2, 2.5] and the choice was made based on expert domain knowledge. The differences between the SLDs can hardly be identified visually, indicating a high applicability of both learning approaches for identifying stable process parameter values as compared to using measured FRFs. However, high deviations between simulated and experimentally determined stability limits can be observed for certain values for the spindle speed. Nevertheless, these deviations could result from inaccurate simulation results and are unlikely a consequence of poor prediction accuracies from the approaches presented in this contribution, since there were no visible deviations between the individual SLDs. The inaccuracy of simulation results can be rooted in several reasons, e.g., neglecting process damping and tool wear or calibrating the compliance models using FRFs, which were measured for a non-rotating tool while the experimentally determined stability limits were a result of analyzing acoustic emission signals of milled slots, obviously using a rotating tool. Moreover, compliance models were only used in the Xand Y-direction. These models were excited using directed force vectors, enabling the calculation of tool deflections based on a weighted interpolation of compliances in Xand Y-direction if the direction of excitation is not exactly aligned with either of these directions. However, the cross FRFs, which represent the relationships between perpendicular directions of excitation and response, were not taken into account. In addition, the solutions of the gradient-based identification of the coefficients of the force model are not unique [74], leading to a limited transferability from calibration experiments to experiments with different engagement situations between the tool and workpiece.

Conclusions
In this paper, an investigation of predicting process dynamics using ML methods for different machine poses is presented. In this context, measurements of FRFs for different poses while using two different machine tools were conducted by utilizing impact hammer tests. An evolutionary-based optimization procedure was used to identify modal parameter values for compliance models that consisted of a set of uncoupled, damped harmonic oscillators and represented the dynamic behavior of the system, consisting of the machine tool, spindle and milling tool, based on an FRF. In order to predict milling dynamics for a given pose, for which no measurements of the frequency response behavior was present, two learning objectives were considered. The first objective was to predict the amplitude and phase for Xand Y-direction, given three pose-dependent features and a frequency value, enabling the capability to predict FRFs. These FRFs were used to calibrate compliance models while using the presented evolutionary-based approach. The second learning objective aimed to predict modal parameter values of compliance models directly, skipping the prediction of FRFs and the subsequently necessary optimization conduction. Both of the learning tasks were evaluated using test sets, including poses with inclination angles that were not used during training. Reasonably low RMSE values could be achieved for all considered learning configurations, except for using linear models for the first learning task. When comparing measured FRFs with predicted FRFs and FRFs which were calculated based on predicted compliance models, predicted FRFs could represent the measured behavior in more detail. However, the prediction of compliance models offers the capability of a real-time evaluation, since the time-consuming optimization task is performed in advance of learning to generate the training set. SLDs were calculated for compliance models that resulted from either measurements, predicted FRFs, or predicted modal parameter values of a set of oscillators using a geometric physically-based simulation system. The comparison between these SLDs showed that there is no notable difference between the different approaches, emphasizing the advantages of predicting the parameter values of compliance models directly for a given pose. The evaluation of simulated SLDs with experimentally determined stability limits showed that simulated stabilities can be used to identify stable configurations to a certain extent. In this context, the accuracy of the used simulation approach limits the applicability of SLDs.
Further research activities could focus on improving the accuracy of the presented approaches for the fusion of peaks of FRFs by incorporating a high amount of observations where this effect is present. Furthermore, the integration of the prediction of compliance models into the simulation system, so that different dynamic properties could be used for each pose along the regarded NC path, could improve the accuracy of simulating different process characteristics, e.g., process forces or tool vibrations. In addition, pose-dependent location errors on the workpiece surface resulting from tool vibrations could be visualized. Moreover, an optimization of process parameter values could be performed while using the prediction of pose-dependent dynamics, to, e.g., identify stable values for the spindle speed for each pose.

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