Probabilistic information fusion to model the pose-dependent dynamics of milling robots

Conventional industrial robots are increasingly used for milling applications of large workpieces due to their workspace and their low investment costs in comparison to conventional machine tools. However, static deflections and dynamic instabilities during the milling process limit the efficiency and productivity of such robot-based milling systems. Since the pose-dependent dynamic properties of the industrial robot structures are notoriously difficult to model analytically, machine learning methods are recently gaining more and more popularity to derive system models from experimental data. In this publication, a modeling concept based on a modern information fusion scheme, fusing simulation and experimental data, is proposed. This approach provides a precise model of the robot’s pose-dependent structural dynamics and is validated for a one-dimensional variation of the robot pose. The results of two information fusion algorithms are compared with a conventional, data-driven approach and indicate a superior model accuracy regarding interpolation and extrapolation of the pose-dependent dynamics. The proposed approach enables decreasing the necessary amount of experimental data needed to assess the vibrational properties of the robot for a desired pose. Additionally, the concept is able to predict the robot dynamics at poses where experimental data is very costly to gather.


Introduction
In order to increase the economic efficiency of milling processes in terms of investment, operating and maintenance costs, conventional industrial robots are an attractive alternative for machining of large workpieces [9,26]. However, the low static and dynamic stiffness of industrial robots often lead to static displacements of the tool or to dynamic instabilities, also called chatter [1,10,21,25]. The static displacements of the tool during the process result in deviations from the target workpiece geometry, whereas dynamic instabilities result in an insufficient surface quality or might even lead to increased tool wear as well as failure of spindle components.
Thus, current research projects address the precise modeling and the system identification of the static and dynamic structural behavior of milling robots. This allows to compensate the estimated errors by choosing a compensated tool path or by choosing stable process parameters [22].
The dynamics of the robot structure are usually formulated in the following, analytic form, depending on the generalized coordinates ⃗ q , the velocities ̇⃗ q and accelerations ⃗ q [14]: M represents the mass matrix, C the Coriolis matrix, N contains gravitational terms and other generalized forces such as joint forces due to the internal joint stiffness and damping properties. are the resulting generalized forces and moments. The efficient simulation of the dynamic behavior of such a rigid body model is provided by software packages, such as the Dynamics Animation and Robotics Toolkit (DART) [12] or the Rigid Body Dynamics Libary (RBDL) [4]. Nonetheless, the identification of the model parameters, such as the inertia, stiffness and damping terms of the joints, remain open issues. Consequently, an incorrect (1) M(⃗ q)⃗ q + C(⃗ q,̇⃗ q)̇⃗ q + N(⃗ q,̇⃗ q) = .

3
parameterization can lead to a potentially inaccurate structural model. In comparison to conventional machine tools, where flexible multi-body simulations are used to model the position dependent dynamics [23], this is especially crucial for robotic milling processes of large workpieces, since the dynamic properties may vary significantly during the process due to the highly nonlinear kinematic structure. Thus, modeling the pose-dependent properties of the robot's structure is still an unsolved problem [2].
Nguyen et al. [15] propose a data-driven approach to model the position dependency of the modal parameters. This approach circumvents the issue of modeling the dynamic properties of a physics-based model accurately, but a large data set is needed to incorporate the full workspace into the data-driven model.
However, a rigid body simulation can still provide a rough estimate of the robot dynamics. Hence, taking knowledge from both information sources into account can alleviate the burden of gathering experimental data in the whole workspace and reduce the necessary amount of experiments. This paper introduces an approach to model and predict the pose-dependent robot dynamics precisely by fusing simulation results with experimental data. The applicability of two information fusion algorithms is demonstrated in comparison to a conventional machine learning approach by inter-and extrapolating the dynamics at the tool center point (TCP) for a variation of the robot's third axis. The concept serves as a feasibility study on how advanced algorithms from the field of machine learning can be fused with physical knowledge to increase the reliability of robot based machining processes.
The paper is structured as follows: the concept is motivated in Sect. 2 on the basis of a comparative study of the vibrational properties determined by a rigid body simulation and the measured dynamic properties. Section 3 provides an overview on methods for information fusion and introduces the main concept of the proposed approach. The model setup and training processes of two information fusion algorithms and a conventional data-driven approach are presented in Sect. 4, followed by a discussion of the results in Sect. 5. Section 6 concludes this work by discussing the present limits and gives an outlook on future research.

Model analysis
In order to visualize the strength of the proposed approach, the dynamics of a KUKA Quantec Prime 240 robot with a high-speed milling spindle were simulated and experimentally measured in dependence of the frequency f ( 5 Hz ≤ f ≤ 30 Hz ) and the robot's third axis z, 3 ( 70 • ≤ z,3 ≤ 120 • ). In both cases, the vibrational behavior is represented by the frequency response function (FRF) H(f , z,3 ) at the spindle close to the TCP, resulting from an excitation in horizontal direction at the same point (see Fig. 1).

Rigid body simulation
The rigid body simulation of the robot dynamics was carried out using Matlab Simscape Multibody with estimated stiffness and damping properties of the joints and estimated mass and inertia properties of the bodies. The model provides three rotational degrees of freedom (DOFs) x , y and z at each joint i with corresponding rotational stiffness and damping properties. Figure 2 illustrates the rigid body model. The model parameters are provided in the appendix.
The frequency response function H(f , z,3 ) is calculated from a linear state space model at a given axis angle of z, 3 . Figure 3 illustrates the frequency response function in the range 70 • ≤ z,3 ≤ 120 • . The simulation was carried out in discrete steps of sim z,3 = 0.1 • , resulting in 501 simulation samples with a frequency resolution of sim f = 0.01 Hz.
As illustrated in Fig. 3, it is clearly visible that the rigid body simulation captured two resonance frequencies: -The first eigenfrequency starts at 7 Hz and increases with increasing axis angle z,3 to 11 Hz.

Experimental data
Similar to the simulation, the FRF H(f , z,3 ) was measured experimentally at discrete axis configurations in the same range 70 • ≤ z,3 ≤ 120 • via impact testing at the driving point illustrated in Fig. 1. The experiment was conducted in discrete steps exp z,3 = 2 • , resulting in 26 measurement samples with a frequency resolution of exp f = 0.0455 Hz (visualized in Fig. 5).
In contrast to the rigid body simulation, the experimentally captured dynamics include three resonance frequencies: -The first eigenfrequency starts at 8 Hz and increases with increasing axis angle z,3 to 10 Hz. The shapes of the first modes of the used milling robot had been experimentally identified in previous works [21,27]. The mode shapes which are assumed to correspond to the three measured modes in this publication are illustrated in Fig. 6. The general pose-dependent behavior of the vibration modes is roughly captured in the rigid body simulation, as the general shape of the pose-dependencies are adequately captured: the pose-dependent behavior of the first mode is well captured in the simulation and the pose-dependent behavior of the second and third measured mode is roughly comparable to the behavior of the second simulated mode. Nonetheless, the simulation does not capture a third mode, but only predicts a second mode with a comparable pose-dependent behavior. It is assumed, that the simulation does not account for the mode which corresponds to the shape in Fig. 6b). The reasons for this issue are twofold: -Inaccurate model parameterization Since the rigid body model is based on a large number of parameters, the identification procedure might have failed to estimate all physical parameters accurately. An inaccurate identification of the mass, stiffness or damping parameters can lead to significantly different dynamic properties. For example, if damping parameters are assumed too high, a mode can be erroneously damped in the simulation. -Unmodeled physical effects Although it is assumed that the measured mode shapes are representing rigid body modes, the simulation might not have taken all physical degrees of freedom into account. Additionally, the simulation does not account for mode coupling, which can be a significant vibrational effect of industrial robots [25].

Probabilistic information fusion
In order to cope with such issues where approximate models can be cheaply evaluated, but precise data is rare or costly to gather, modern machine learning algorithms are capable of fusing information from different data sets and information sources with different fidelity levels.
Hence, such algorithms are also referred to as multi-fidelity information fusion algorithms [5,16].
Meng et al. [13] proposed a deep learning approach, where two artificial neural networks are coupled hierarchically to fuse information from two data sets with a different fidelity level. An uncertainty estimation is included based on a Dropout approach [7]. Additionally, since the implementation of this approach makes use of automatic differentiation methods, the training objective can also take physically motivated bounds in an analytic form into account, as previously published in [19,20].
Similarly, the information fusion can be based on probabilistic, Bayesian inference [3,6,11,17,18]. The approaches mainly make use of Gaussian and deep Gaussian process regression techniques and differ in the ability to fuse data with linear or nonlinear correlations. In contrast to deep learning-based approaches, Bayesian inference algorithms rely on a statistically well-founded theory. Thus, they provide a reliable uncertainty estimation.
Such a probabilistic information fusion scheme can be set up to infer the mapping between a simulation of the robot's structural dynamics and the experimental data: The pose-dependent dynamic behavior can be easily approximated by a rigid body simulation. In contrast, experimental data, for example gathered via impact testing or automated shaker experiments, is costly to generate, but provides a precise representation of the pose-dependent dynamics. The approach presented in this paper addresses this issue by fusing the information of both sources using a multi-fidelity information fusion algorithm.
The setup of such a multi-fidelity information fusion scheme is described more formally in the following: As previously introduced, the general objective is to find an unknown mapping between an input space X and an output space Y . The latter is assumed to be one-dimensional: To find such a mapping, a training data set D is given, which consists of m input data samples ⃗ x ∈ X with corresponding output data y ∈ Y: In the present case, the objective is to find a mapping between the input data ⃗ x = [f z,3 ] T and the output data space y = H . For the depicted scenario, there are two information sources, which provide information on this unknown relationship. The goal of the information fusion algorithm is to find a suitable mapping between the approximative but cheap data D LF (also called low fidelity data) and the precise but costly data D HF (also called high fidelity data).
If there is a linear relationship between the two fidelity levels, the relationship between low fidelity LF (⃗ x) and high fidelity HF (⃗ x) can be described with a constant scaling factor c and a bias term (⃗ x) [3]: In case of a nonlinear relationship between the two fidelity levels, a constant scaling factor is not sufficient. Instead, the functional relationship between low and high fidelity data must incorporate a nonlinear transformation nl ( LF (⃗ x), ⃗ x) [3]: For the given scenario, the two nonlinear information fusion algorithms NARGP [17] and MFDGP [3] are used to infer the possibly nonlinear mapping between a simulation of the robot's structural dynamics and the experimental data. Both approaches have a hierarchical, layer wise (deep) architecture of Gaussian processes in common, where a separate Gaussian process is conditioned for each fidelity level. The approaches differ in terms of which data set is used to condition each Gaussian process and in their specific training procedure.
In the following, the setup and the training processes of both information fusion algorithms NARGP and MFDGP are described. For comparison, the setup and the training process of a conventional, Gaussian process model are also described. In Sect. 5, the performance of the information fusion approaches is assessed. The conventional Gaussian process regression model serves as a benchmark to illustrate the superior performance of the information fusion techniques.

Model setup and training
To illustrate the performance of such a data-driven modeling approach, the data set D HF is divided into a training data set D HF,train and a test data set D HF,test . The training data set D HF,train consists of the first 70 % of the data regarding the axis angle z,3 , which corresponds approximately to the range of 70 • ≤ z,3 ≤ 105 • , whereas the test data set D HF,test consists of all remaining data samples in the range 105 • < z,3 ≤ 120 • (illustrated in Fig. 7). D LF,train consists of 12801 data points, while D HF,train only consists of 2622 data points. Both data sets are taken equally distributed from the data shown in Sect. 2.1.
In the following, it will be shown, that both information fusion algorithms are able to combine the information from both data sets D LF and D HF,train . The model setup and training process was conducted with the Python library Emukit [24].
As proposed by Perdikaris et al. [17], the most general NARGP model is based on two consecutive Gaussian processes with a radial basis function kernel k(⃗ x, ⃗ x � ): where k(⃗ x, ⃗ x � ) is the kernel to model the covariance between two data points ⃗ x and ⃗ x ′ , 2 is the variance and l the length scale of the kernel.
The MFDGP model incorporates a more complex kernel design. This layer-wise kernel design is described in [3]. The radial basis function kernel was extended in such a way, that linear relationships between the low and high fidelity data are better addressed. Although different kernel designs are possible, the model setup in this work follows the proposed kernel design of Perdikaris et al. [17] and Cutajar et al. [3], since physical expert knowledge is deliberately not included by manual kernel shaping.
To speed up the training process, MFDGP relies on a sparse variational approximation method with 800 inducing points and a mini-batch size of 50. Comparable to the approach of Cutajar et al. [3], the model training was conducted consecutively: in the first step, the model was trained with a fixed variance, followed by a second training step which includes the model variance in the training process. Table 1 summarizes the model and training parameters for the MFDGP model.
As previously mentioned, a conventional Gaussian process regression model can serve as a reference benchmark. The model training was conducted using the Python library Gpy [8]. The Gaussian process regression model was set up with an additive standard radial basis function kernel and a bias kernel to shift the mean as follows: In the following section, the training results are compared and the model performance of each algorithm is assessed.

Results
The proposed information fusion approach can improve the model accuracy significantly in comparison to conventional machine learning approaches. Figure 8 shows the prediction results for the training and the test data set. Since all three algorithms provide probabilistic models, the frequency response function H(f , z,3 ) is not only represented by the expected value (f , z,3 ) , but also by an uncertainty estimation based on the standard deviation (f , z,3 ).
The model performance can be evaluated by considering the model's capabilities to interpolate and extrapolate the dynamic properties of the robot: It can be observed, that the conventional Gaussian process model is capable to represent the robot dynamics for axis angles where experimental data is nearby. As seen in Fig. 8a), the model is able to interpolate within the training data set (axis angles between 70 • and 105 • ), but the model's accuracy quickly deteriorates with increasing axis angle in the test data set (from axis angle 105 • onwards). Thus, a conventional data-driven approach is unable to extrapolate the robot dynamics, which results in poor predictions of the  (f , z,3 ) , as the expected value GP (f , z,3 ) tends to a mean value (see Fig. 8d). The MFDGP model can extrapolate the dynamics better that the conventional Gaussian process regression model (see Fig. 8b). Nonetheless, similar issues arise: the MFDGP model lacks a precise extrapolation of the frequency response, while the uncertainty increases with increasing axis angle in the test data set from axis angle 105 • onwards (see Fig. 8e).  In contrast, the trained NARGP model is able to extrapolate the pose-dependent dynamic behavior (see Fig. 8c). Since the amplitudes of the measured frequency response functions vary especially at the resonance frequencies due to a nonlinear structural dynamic behavior, the NARGP model is able to incorporate this knowledge implicitly in its uncertainty NARGP (f , z,3 ) (see Fig. 8f).
In the following, the model accuracy is assessed in detail. Figure 9 illustrates the accuracy of a NARGP model for the training as well as the test data. Additionally, the simulation accuracy is illustrated as a benchmark. As indicated by the black dots, the NARGP model performs very well on the training data set. Additionally, it is visible that the model also performs well on the test data set (indicated by the blue dots), whereas the rigid body simulation accuracy is significantly worse, since it does not incorporate all vibrational modes. Furthermore, the NARGP model's standard deviation provides a comprehensible uncertainty quantification. Table 2 summarizes the model performance by quantifying the prediction accuracy of the rigid body simulation and the three models, pointing out the significantly more accurate prediction of the NARGP information fusion scheme in regard to the coefficient of determination R 2 and the rootmean-square error RSME.
In order to assess the interpolation capabilities and the efficiency of the depicted approach using the NARGP algorithm, the model performance is evaluated by gradually reducing the number of experiments used in the training data set. For each assessment i ∈ {1, 2, 3, 4} , only every ith experiment of the whole data set was used for training as high fidelity data. On the contrary, the complete simulation data was used as the low fidelity data set for each training.
The results, illustrated in Fig. 10, indicate that the prediction accuracy decreases with fewer experiments used for training, especially regarding the second and third measured eigenfrequency. Nonetheless, the prediction remains highly accurate, even when the number of considered experiments is reduced from 26 to 9 measured frequency response functions (i.e. a reduction of over 65%).

Conclusion
In this paper, an approach to improve the modeling accuracy of structural dynamics models of milling robots was presented. This approach is based on an information fusion algorithm to combine physics-based rigid body simulation models with experimental data. The results show, that it outperforms both, the rigid body simulation and conventional machine learning approaches in estimating and extrapolating the pose-dependent robot dynamics. It is worth mentioning that the presented approach is not only limited to model the pose-dependent of dynamic properties of industrial robots such as milling robots, but could also be used to model the position-dependent dynamics of conventional machine tools.
Nonetheless, there are open issues which need to be addressed in future research: -Economic utilization of the prediction uncertainty An essential advantage of Bayesian approaches, such as the used deep Gaussian processes, is the well-founded incorporation of a statistically motivated model uncertainty.
It remains an open issue how the prediction uncertainty can

Fig. 9
Prediction accuracy assessment for the NARGP model in comparison to the rigid body simulation; the prediction of the NARGP model is marked in black and blue for the training and test data sets. The uncertainty of the NARGP model for the test data set is highlighted in blue be incorporated into the robust design and parametrization of robotic milling processes. -Computational complexity As of now, the power of this approach has only been illustrated regarding a onedimensional robot work space (as a variation of the third axis). In order to make use of such models in industrial applications, the depicted approach needs to be extended to more degrees of freedom, which results in significantly larger data sets. Nonetheless, the computational complexity of the given algorithm scales nonlinearly with the number of training samples (see [3]). Thus, it is necessary to examine other feature space representations. In addition, reducing the computational complexity also enables more sophisticated validation schemes, such as cross validation, which is favorable when the approach is extensively used in larger workspace areas. -Nonlinear simulation model The depicted approach has been illustrated using a linear, pose-dependent state space model of the robot's structural dynamics. Since the information fusion heavily relies on an accurate estimation of the most important vibrational properties, the physics-based model should also incorporate nonlinear effects.
Acknowledgements Open Access funding provided by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.