Optimal design for compliance modeling of industrial robots with bayesian inference of stiffnesses

In this paper a cost and time efficient approach to setup a compliance model for industrial robots is presented. The compliance model is distinctly determined by the gear’s stiffness parameters which are tuned by an optimal design of experiments approach. The experimental setup consists of different poses of the robot’s axes together with the applied force at the tool center point (TCP). These robot poses represent together with defined forces the experimental setup where the deviation of the robot under defined force is measured. Based on measurements of the displacement of the TCP the stiffness parameters for the compliance model are estimated and afterwards validated in new experiments. The efficiency of this approach lies in the reduced amount of experiments that are needed to identify the stiffness parameters that are parameters inherent to the compliance and the less complex experimental setup.


Introduction
Industrial robots have major advantages over machine tools in terms of flexibility and workspace. The flexibility of industrial robots is expressed both in the accessibility of orientations of the TCP in the workspace as well as the ability to perform a wide variety of manufacturing processes. Also the capital costs for industrial robots are up to three times lower than the costs of a machine tool [1]. Therefore industrial robots are more and more treated as an alternative to machine tools.
However industrial robots have lower working accuracy compared to machine tools because they are less rigid [2]. Compared to a machine tool the absolute accuracy of an industrial robot can be between 20 to 50 times lower especially under changing process forces [1]. The lower accuracy results mainly from geometric erros and a low stiffness. Therefore industrial robots are suitable for tasks where the expected forces and the expected accuracy are low [3]. To compensate or minimize the main geometric errors that result from manufacturing and assembly of the industrial robot a producer specific calibration is carried out and deposited in the control of the robot before the robot is shipped [4]. Therefore, after compensating the geometric erros, the low stiffness of the robot is the main reason of the low accuracy. The research on stiffness improvement mainly concentrates on three aspects: stiffness compensation, posture optimization and robot structure optimization [5]. The low static stiffness, that has been identified in several works [6,7] as the significant error influencing the working accuracy of the industrial robot especially under process forces, is mainly caused by the used gears and their low rigidity. The influence of the low rigidity of the gears outweighs the deformation of the structure and the tilting rigidity of the bearings [6]. This results in a pose dependence of the stiffness of the industrial robot [7]. There are different approaches to compensate this stiffness behavior such as constructive changes in the robotic structure that comes with reduced mobility and flexibility of the robot. Another approach is to model and predict the compliance behavior of the industrial robot under measured or predicted process forces and compensate the deviation with offline or online compensation [1].
The effects of the compliance behavior of the industrial robot especially under process forces has been examined and modeled for example by [1,7,8]. The models are based on the forward transformation that is extended with the elastic deformation of the joints [8]. In this modeling approach the elastic structures are modeled by rigid bodies connected by virtual elastic joints. These virtual joints are then implemented in the forward transformation [1]. The level of complexity of these models differs depending on whether only the torsional rigidity is considerated or also the tilting rigidity of the structure [1].
Since the deviation of the TCP mainly results from the limited torsional stiffness of the gears often only the torsional rigidity is implemented and the links are considered to be almost stiff [4]. For the forward transformation the kinematic description of the robot is necessary, to describe the posedependent compliance of the robot the stiffness parameters of the gears are necessary. To define the robot specific stiffness parameters, mostly complex setups were used for example in [1] the stiffness parameters were determined by applying observable cubes to the robot structure. By determining the displacement of the cube while a defined test force was applied on the TCP the compliance was calculated. To apply the force a cable pull system was used and, to measure the determination of the parts of the robot, a laser tracker was used. This example shows that the determination of the stiffness parameters is time and resource consuming.
Cordes [6] used a less resource consuming approach by measuring the global compliance of the robot structure. In twelve different positions a defined force was applied on the TCP and the shift of the TCP was measured by cost intensive metrology. By combining the measurements and minimizing the sum of squared residuals of the over-determined system the parameters were determined.
In this paper a more cost and time efficient approach is implemented which is based on optimal design of experiments for variance-minimal Bayesian inference of the gear stiffness parameters. The experimental setup consists of the angular poses of the robot's axes and the direction as well as the magnitude of the force which is applied on the TCP. In the experiments the quantity of interest is the deviation of the TCP which is measured by laser triangulation sensors. Since the measurement process is subject to aleatoric and epistemic uncertainty [9,10] so are the estimated parameters. The experimental setup is optimized so that the uncertainty of the estimated stiffness parameters which is measured by a scalar function of the covariance matrix becomes minimal. In the second step measurement data is collected within the optimized experiment and infer the compliance parameters in a Bayesian setting [11]. The prior is chosen from previous estimations in [6] and enters the estimation process.
Optimal design of experiments is a broad field of research and describes an experimental design that is optimal regarding statistical criterion [12][13][14][15][16]. In general, the literature can be classified into topics concerning optimal sensor placement [15,[17][18][19] and optimal choice of inputs [16,[20][21][22][23][24]. In this case the sensors which observe the movement of the TCP are fixed but the angular poses of the axes and the applied force can be varied. Then it is refered to these conditions as the inputs from now on. This paper is structured as follows. In Sect. 2 the compliance model is described in full detail. In Sect. 3 the approach to optimal design of experiments for Bayesian inference of the gear compliance parameters is presented. The experimental setup is presented in Sect. 4 and the numerical results are shown in Sect. 5 before a conclusion is given in Sect. 6.

Compliance model
The aim of the compliance model is to calculate the deviation of the industrial robot caused by process forces [6]. To describe the translational deviation in the cartesian space the elastastic behavior of the gears are combined with the forward transformation. The forward tranformation is the foundation of the model and describes the position and orientation of the TCP in the cartesian space based on the joint angles i [8]. To set up the model, on each joint of the robot a coordinate system is equipped as shown in Fig. 1. The relation between the coordinate systems and therefore the robot parts are described by tranformation matrices according to the Denavit Hartenberg convention [25]. For the complicance model this tranformation has to be changed so that the transformation ends with the rotation around the Z axis [25]. The transformation from one coordinate system to another is based on rotation and translation parameters.The rotation parameters are i which describes the rotation around the X axis and ̃i which describes the rotation around the Z axis. The translation prameters are a i which describes the translation in X translation and d i for the Z translation [1] which yield the forward transformation (1) The rotation ̃i includes the joint angle i [25].
In this approach it is assumed that the links of the robot are rigid and the joints are elastic in direction of rotation with high radial stiffness behavior. Therefore the deformation of the gears is modeled with the behavior of linear elastic torsional springs [4]. To describe the deviation of the robot the linear elastic behaviour of the gears is used to extend the forward transformation by an additional transformation that is based on the joint stiffness as in [8] described.
The rotational part of the matrix describes the deformation of the virtual joint due to a measured force at the TCP with the angle Δ that describes the torsion around the Z axis of the gear due to the force at the TCP [1]. To calculate Δ the equation for linear torsion springs is used where M i represents the torque that is caused by the cutting forces and k i the stiffness parameters of the joints. The torques can be calculated by building the cross product of the vector processing forces ⃗ F and the transformation matrix towards the specific coordinates systems of the axis ⃗ T i . For this the forces engaging at the TCP ⃗ F TCP have to be transformed into the joint coordinate system with so that For the stiffness parameters k i it is assumed that the joint behavior is linear and therefore the stiffness parameters are assumed to be fixed parameters [6]. With T R (Δ i ) an additional tranformation matrix with rotation around Z axis is implemented in the forward transformation. The conventional forward kinematics provide the unloaded TCP position ⃗ TCP ideal , the extended forward transformation the TCP position under process forces ⃗ TCP F [25]. The force-induced deviation can be calculated with The advantage of this analytic description of the deviation is the fast calculation time compared to simulations [25].
For this work the six-axis heavy duty robot ABB IRB 6660-205/1.9 is used. It is build for milling processes with a closed-loop mechanism in form of a parallelogram which can be seen in Fig. 1. This causes a coupling between the joint angles 2 and 3 which is described in the model as To implement the model for this robot the parameters i , a i , d i and k i have to be determinded. While the parameter i , a i , d i can be extracted from the documentation of the robot or through measurements of the robot structure [26], the stiffness parameters have to be determined in experiments. To determine the stiffness parameters an optimal design of experiments approach for Bayesian inference is used.

Bayesian inference and optimal design of experiments
In optimal design of experiments the best experimental setup is chosen such that the desired parameters in this setting the stiffness parameters can be estimated from data with minimal variance. In this setting, the orientations of the robot joint angles and the direction and magnitude of the force f are optimized such that the uncertainty in the estimated stiffness parameters k becomes minimal. The input variables u ∶= ( , f ) are introduced and used from now on. First the Bayesian framework is described and then the optimal design of experiments approach. The compliance model ( k , u , ) maps the parameters k , the inputs u and some constants to the three-dimensional deviation of the TCP. This is the quantity of interest in the experiments which are measured via lasers. The compliance model serves as a computer model which needs to be calibrated and validated. From a probabilistic point of view the data z is assumed to be a realization of a random variable Z with a specified distribution. It is assumed here that the data is subject to Gaussian observational noise: where z ⋆ is the true but unknown quantity of interest. If this model is correct than it describes the data well for all inputs u from an admissible domain U ad : where k ⋆ are the true but unknown stiffness parameters. The aim is to find a good approximation k of k ⋆ from a Bayesian inference approach. For a given realization z of the data Z and for a given realization of the noise E it is required that Within a Bayesian approach, the estimated parameters have a posterior probability distribution ( k ‖ z ) that is derived from the data likelihood ( z ‖ k ) and the prior probability distribution 0 via the Bayes formula. Let 0 ∈ N ( k 0 , 0 ∕ ) be a Gaussian prior, where > 0 .
Moreover, let (⋅) be the density function of the noise random variable E . In view of Eq. (9), the data likelihood ( z ‖ k ) has the density function z ↦ ( z − ( k , u , )) . Then the posterior ( k ‖ z ) is given by where C is a normalization constant and ‖x‖ B ∶= √ x ⊤ Bx is the weighted norm of a vector x with a compatible matrix B . The maximum a posteriori estimator (MAP) is defined to be a point k that maximizes the posterior: compare [11]. It is evident that the estimated parameters depend on the inputs u and the data z , i.e., k = k ( u , z ).
The compliance model is nonlinear in k . Thus, the posterior is not Gaussian. However, the model is linearized at the MAP point k to obtain Gaussian posterior probability distributions. In the linearized case the posterior covariance matrix is computed [15]: where J( k , u ) is the sensitivity matrix of the model with respect to the parameters: Let 1 − , where ∈ (0, 1) , be a confidence level. Then the confidence region G( k , C post , ) around the MAP point k is an n k -dimensional ellipsoid: (1 − ) is the quantile of the chi-squared distribution with n k degrees of freedom, see [12,16].
Optimal design of experiments aims at minimizing the posterior covariance matrix C post . In order to avoid illposed optimization problems [12], one minimizes a scalar design function Ψ(⋅) of the covariance instead. The following design criteria are frequently used [13]: In this paper, Ψ(⋅) = Ψ A (⋅) is chosen, i.e., the trace of the posterior covariance matrix is computed. In order to reduce the complexity of the problem and since the measurement data z is not available at this point k ( u , z ) ≈ k 0 is fixed. Thus, to minimize the trace of the covariance matrix The following constrained optimization problem arrises: In our setting the set of admissible inputs U ad is determined in the following way: The first condition states that the joint angles and the components of the force vector have to lie within a minimum u min and a maximum u max range. The ranges are defined in Sect. 4. The second condition states that the orientation of the TCP is within a range that can be reached by the robot and is suitable for milling processes. Note, that the functions g and h are nonlinear. Thus the problem (18) is restated more explicitly: This optimization problem is solved with a standard Sequential Quadratic Programming (SQP) solver with BFGS updates, see [27].

Experimental set up for stiffness parameter identification
Based on the compliance model which was developed in Sect. 2 and the in Sect. 3 introduced aproach the experimental setup for the stiffness parameter identification has been implemented. The experiments are carried out on the six-axis heavy duty robot with a spindle as an effector. The robot cell that includes the heavy duty robot is built to carry out the hybrid manufacturing chain including additive and subtractive processes. Also included into the cell is the force sensor Omega 160 IP65 from ATI that is placed between the process effector and the sixth joint. The robot cell also includes a work table, an effector docking system to change the spindle for the additive manufacturing head and measurement equiment to define the position of the tool and machining parts. The robot cell is shown in Fig. 2.
The setup to identify the deviation of the TCP under force contains the application of a defined force and the measurement of the deviation. To apply the force a pneumatic cylinder is used. It contains a force sensor to measure the applied force on the TCP so that the applied force can be adjusted. As in Sect. 3 described the force vectors are represented in the first condition of the inputs. The input u min is defined as −300 N and the u max as +300 N due to the experimental setup. To measure the deviation of the TCP under the defined force, three laser triangulation sensors of the type Keyence LK-G32 are used. The sensors have a measuring range ±5 mm and an accuracy of 1 μm . They are applied around the TCP and measure the three directions X,Y and Z.
To measure the force that is applied in the robot coordinate system also the integrated force sensor of the robot is used. With this the vectors of the applied force can be defined in the robot system. The experimental setup can be seen in Fig. 3.
The advantages of this experimental setup are that it is a simple and fast setup, no changes of the robot system are necessary and the experiment does not influence any existing  compensation or calibration. The measurements are executed for one position that is elected considering the workspace the robot operates in. To be able to represent the static force conditions of the peripheral milling processes, the experiments are examined with forces applied in X and Y direction. For the measurements only the orientation of the effector has to be changed and the direction of the applied force, the position in the cartesian coordinate system is not changed. The orientations of the effector and the applied force are calculated with the in Sect. 3 developed model. One measurement consists of one orientation of the effector and an applied force in one direction. To ensure that the orientation of the effector calculated in the model is reachable, the solution space of the model is limited. The orientation of the TCP is only allowed to be within the angle of ±80 • between the direction of the vector of the joint six to the TCP and the direction vector between the robot base and the TCP regarding the Z plane. is here defined as: The in Sect. 3 developed model needs initial values for the stiffness parameters to calculate the optimal joint angles i . Therefore the stiffness parameters from [6] are used as can be seen in Table 1.
The calculated design of experiments contains one position in the cartesian coordinate system, four orientations of the effector to reach this position and a total of five measurements. Three measurements are to be carried out with applied force in Y direction and two in X direction of the robot base coordinate system. The data of all sensors are collected via matlab and synchronised. Afterwards they serve as input variables for the model created in Sect. 3.

Parameter identification and validation
For the identification of the stiffness parameters different sets of data were used to analyse the influence of the number of experiments. The data of one measurement in one direction, of three measurements in the same direction and the data of all measurments were used to calculate differents sets of stiffness parameters. The results were then validated through comparison with the measured deviation in validation measurements. Also the set of stiffness parameters from the literature shown in Table 1 was compared to the calculated sets of stiffness parameters. An extract of the calculated sets of stiffness parameters from different measurements are shown in Table 2.
It can be seen that some sets of stiffness parameters differ enormus from the initial parameters from [6] shown in Table 2 while others correspond with the initial parameters. This can be explained by the direction of the applied force. For example during the measurements in the X direction the joint six was not loaded so that no information about this joint is contained in the data of this measurement. To validate the parameters and indentify the best set of stiffness parameters the deviation of the TCP in various positions and orientations of the main workspace of the robot was measured. Figure 4 shows the measured deviation and the with different stiffness parameter sets calculated deviation for one validation position. The results of the other carried out validation experiments for different positions are similar to the results shown in Fig. 4. It can be seen that overall the deviation of different positions and orientations can be predicted. The accuracy regarding the deviation in Y direction is higher than in X direction which can be explained by the facts that in Y direction more measurements were executed to identify the set of stiffness parameters, the deviation is higher and more gears were loaded during the experiments. The figure also shows that the results of the calculated deviation based on the calculated stiffness parameters are closer to the real deviation than the calculated deviation based on the set of stiffness parameters from literature. It can be seen that the stiffness parameters from measurements in Y direction and the stiffness parameters from all measurements have the highest accuracy while the stiffness parameters from the literature lead to the lowest accuracy. Overall the parameters from measurements in Y direction match the most if the outliners are disregarded. This can also be seen in the covariant matrices that were calculated together with the stiffness parameters see Sect. 3. The covariant matrix of the measurements with the force direction X have the highest values while the covariant matrices of the measurements in force direction Y and of all measurements are close. They are by the factor 10 lower than the values in the covariant matrix of the measurements with the force direction X. The covariant matrix of the measurements in force direction Y are a little lower than the one of all measurements. The results of the covariant matrices are reflected in the results of the calculated deviations in Fig. 4 compared with the measured deviations. The stiffness parameter sets with the lowest covariant matrices results calculated are closest to the measured deviations.
Overall the data show that the measurements produce stiffness parameters with which the deviation can be estimated more accurate than with the stiffness parameters from the literature of the same robot model. Also the results show that the in Sect. 3 implemented approach is valid and with the calculated covariant matrices an estimation of the quality of the results can be made. Compared to [6] the amount of measurements are reduced up to 75 % and the used experimental set up is less cost intensive and les complex compared to [1,6]. Therefore this approach helps to identify the stiffness parameters by using a simple and cost efficient experimental setup.

Conclusion and outlook
The deviation of the TCP under process forces is the main influence on the working accuracy of the robotbased milling process. It limits the fields of applications of these systems. To extend the application possibilities the deviation of the TCP under static process force has to be detected and predicted. This paper shows an aproach to predict the deviation of an industrial robot under static process forces by identifying the stiffness parameters of the joints based on the bayes approach. Based on this approach an optimal design of experiments for the measurement of the deviation under a specific force is calculated. With this approach the uncertainties of the results are reduced and the level of information a measurment contains is increased. This approach leads to less measurements of the deviation for the identification of the stiffness parameters of an industrial robot which is more time efficient as other approaches. Also the experiments itself are less complex due to a simple setup and therefore more cost efficient than other approaches. The results show that the identified stiffness parameters lead to a sufficient accurate prediction of the deviation and thus the identified stiffness parameters are more accurate than the stiffness parameters from the literatur for the same robot model. The advantages of this experimental setup are that it is a simple and fast setup, no changes of the robot system were necessary and the experiment did not influence any existing compensation or calibration. Also the quality of the The experimental setup only allowed the application of force in one direction during a measurement. To increase the accuracy of this approach it is assumed that in an ideal experimental setup forces in two or three directions can be applied and with that the amount of measurements will be further reduced. Therefore the accuracy of the estimated set of stiffness parameters could be increased.
Nevertheless the results show that even with a less complex design of experiments the stiffness parameters are sufficient for the use for online or offline compensation. With the deviation model and the identified stiffness parameters an offline compensation will be implemented. Based on the forces of previous milling processes the deviation will be predicted and the new trajectory can be calculated. The paper shows that with this approach only a simple design of experiments is needed to identify sufficient accurate stiffness parameters to predict the deviation of the TCP. With this approach the influence of the process forces on the deviation during the milling process can be predicted and in a second step reduced so that the quality of the produced parts will be increased.