Error Prediction in Industrial Robot Machining: Optimization Based on Stiffness and Accuracy Limit

Among the advantages of using industrial robots for machining applications instead of machine tools are flexibility, cost effectiveness, and versatility. Due to the kinematics of the articulated robot, the system behaviour is quite different compared with machine tools. Two major questions arise in imple-menting robots in machining tasks: one is the robot’s stiffness, and the second is the achievable machined part accuracy, which varies mainly due to the huge variety of robot models. This paper proposes error prediction model in the application of industrial robot for machining tasks, based on stiffness and accuracy limits. The research work includes experimental and theoretical parts. Advanced machining and inspection tools were applied, as well as a theoretical model of the robot structure and stiffness based on the form-shaping function approach. The robot machining performances, from the workpiece accuracy point of view were predicted.


Introduction
Among the benefits in applying robots to machining tasks that were first reported in the 1990s were increased flexibility and lower costs. Industrial robots (IRs) integrated in cutting applications are primarily used for prototyping, cleaning, and pre-machining of cast parts as well as end-machining of middle tolerance parts [1]. Notwithstanding their advantages, IRs suffer from inherent weaknesses under the presence of machining process forces, namely: low positioning accuracy, vibration, and deflections. For articulated robots, repeatability is inherently dependent on its reach distance, such that the greater the reach distance, the lower the repeatability [2]. Therefore, in robotic milling applications, the process final results are unwanted trajectory deviations, which lead to errors in target dimensions and reduced surface quality of the workpiece (WP). These deviations are mainly caused by static offset overlaid with low frequency tool oscillation [3]. Today's optimum repeatability levels for industrial articulated robots can reach ±0.03 mm [4], which is sufficient for many low-to medium-accuracy part machining tasks. Robot positioning accuracy and repeatability were investigated and reported in [2]. The tool center point (TCP) position error result for the KUKA KR 16 -2 robot without loading during movement along a cube with a side length of 1000 mm was 0.15 mm. Previously reported TCP deviation from the nominal programmed position using the LOLA 50 robot under a static load of F ystat = 100 N was 0.3 mm [5]. Additional deflections of 0.25 mm were measured under milling loads of 100 N using the KUKA KR 210 robot, which was consistent with the expected compliance [4]. The accuracy of robots in machining was reported for various operation types and under several processing conditions, for example, a deviation error of 0.19 -0.55 mm was obtained for the KUKA KR 125 robot in milling with a 300 N load [6]. In another milling test using the IRB 6640 robot with a 500 N load, an error of 1 mm before compensation and 0.4 mm after compensation was achieved [7]. An error of 0.2 -0.35 mm was reported in a drilling test using the KUKA KR500-2 robot with a process load of 1000 N [8]. Although studies on the positioning errors, measurement, and compensation, in robotic applications, were conducted since the 1990's [9] [10], it has been difficult to improve the accuracy of machining tasks. The main reason was the lack of a reliable programming tool for predicting possible TCP position error in relation to the WP geometry and to the robot workspace.
In this paper, the developed method for enabling the end-user to predict, through simulation, the achievable accuracy in robot-aided machining and to optimize the WP position in the workspace in terms of the achievable machining accuracy were described. The experimental validation procedure was performed. The developed software, based on the form-shaping function (FSF) approach and previously applied for modelling machine tools [11] [12], is used here for the first time for articulated robot research.

Overview of Research Approach
A flowchart illustrating the research approach, which comprises theoretical and experimental stages, is shown in Figure 1. The form-shaping function (FSF) approach made it possible to construct the kinematic model of the investigated serial robot. (Steps A.1, A.2). By solving the inverse kinematic equations, the robot spatial configuration, including six values of the joint rotations were defined (Step A.3). A set of joint compliance measurements were then applied Engineering (Step B.1). Based on the joint stiffness/compliance model, the joint compliance matrices were composed (Step A.4). Then the full compliance matrix of the robot was calculated (Step A.5). A two-step procedure using the Renishaw Q10 ball-bar system was applied in order to teach the robot TCP the planned movement along a circular path (Step B.2). A ball-bar test, which was run while the predefined circular TCP path was being performed, determined the positioning accuracy and repeatability of the robot when no external forces were applied (Step B.3). Experiments, including circular groove milling, were performed. A high-speed Mini-spindle, and the Kistler dynamometer system were used for machining and measuring the cutting forces during cutting (Step B.4). The accuracy of the machined circular groove, roundness and position was measured using of-line CMM-XYZ (Step B.5). In order to calculate the tool tip deviations from the nominal path during cutting, the measured cutting forces were used as external data in the simulation model (Step A.6). The simulation model was verified by comparing the results of the deviations measured in the experimental part, i.e., deviations from the nominal circular groove, and the results obtained with the simulation software (step C).
Simulations were run to obtain a visual description of the TCP deviation vectors, to map the deviations, and to define the workspace within a predefined deviation limit (Steps A.7, A.8, A.9). In Step D, concluding this study, the achievable workspace within a predefined accuracy limit is given. Figure 2 depicts the proposed model of 6R robot, which considered as serial kinematic chain with six kinematic links, revolute joints, labeled from base to tool, as follows: S, L, U, R, B, and T. The positive direction of rotation for each joint is also shown in Figure 2.

Form-Shaping Function of the Serial Robots
The applied modeling was based on the error budget theory and FSF approach [11] [12]. The results of this study enabled us as follows:  To predict the achievable tolerances on a machined workpiece (WP).  To optimize WP mounting position relative to robot configuration.  To define the permissible mechanical loadings, and, therefore, the permissible machining parameters (tool geometry, feed rate, spindle speed, and depth of cut).
The form-shaping system (FSS) of the machine tool or robot consists of an ordered aggregate of machine links, whose relative positions and mutual movements ensure the specified travel trajectory of TCP with respect to a WP.
The main mathematical model of the FSS theory is the FSF. Depending on the object to be transformed, a position vector, an orientation vector, or a screw, three types of the FSF are considered.  x y z = r , In Equation (1), x 0 , y 0 , and z 0 are the coordinates of an FP referring to the frame S 0 whereas, x n , y n , and z n are those referring to the frame S n ; and 0 A n is the 4 × 4 manipulating matrix of the FSF presenting a product of the cofactor-matrices

Kinematic Model of 6R Robot
The manipulation matrix 0 A 15 is the product of the fifteen 4 × 4 matrices given in the last column of Table 1 [13].
where descriptions of all cofactors included in the expression for the manipulation matrix (Equation (6)) are shown in Table 1. The cofactor-matrices in Equation (6) presented in Table 2 [12]. The kinematic model of the 6R robot associated with Table 1 shown in Figure  3.

Solving the Inverse Kinematic Problem
The inverse kinematic problem is to find solutions of equations for the FSF system (Equation (4)), which consists of six equations with six unknowns, one for the rotation angle of each of the six joints: Engineering φ: rotation angle of joint S relative to the base around the Z-axis of the base coordinate system.

Stiffness/Compliance Model of the Joint
The robot links assumed as rigid bodies, and the joint stiffness (which includes control loop stiffness and actuator mechanical stiffness) represented by a linear torsion spring. The deformations of joint bearings that theoretically presented as springs with nonzero compliance were not considered in this study.
The compliance values of the joints (which include the compliance of the control loop and the mechanical compliance of the drives) were measured previously for the 6R robot type MH-12 YASKAWA by loading each link with a known external load and measuring the angular deviation of the current link relative to the previous one. Since the presented study has the estimated character of the achievable accuracy of robot machining, the compliance of the joints was measured approximately, and the process of these measurements is not de-

Stiffness/Compliance Model of the Serial Robot: Full Compliance Matrix
As shown in [12], the 6 × 6 stiffness-compliance matrix of the terminal link of a serial robot with respect to the base, expressed through 6 × 6 joint stiffness matrices, is as follows: where C 0N and K 0N are the full compliance (stiffness) matrix of robot, K ink is the 6 × 6 stiffness matrix of the kth link relative to its coordinate system, and C ink is the 6 × 6 compliance matrix of the same link. All components of Equation (7) must be represented in the same coordinate system. To obtain the full compliance matrix C 0N by using Equation (7), we need to calculate all compliance matrices of all joints relative to the platform coordinate system using Equation (8): where the definition of C ink is given in Table 3. The Equation (7) and Equation  T T  T  T  T  T  T  T   T  T  T  T  T T T  T  T  T  T  T  T  T   T  T  T  T  T  T  (9) where one of six one-parametric matrices T j (q k ), 1, 2, , 6 j =  , describing a coordinate transformation of a screw from the frame S k associated with the kth link of the chain into the frame S k−1 , is defined in Table 2. Equation (7) for calculating the full compliance matrix is: where each joint compliance matrix calculated by substituting the corresponding Equation (9) and expression from Table 3 into Equation (8).

Calculation of the Robot Tool Deviation under Forces Acting during Milling Process
The forces F x , F y , and F z were measured during milling operation of a circular furrow on aluminum plate, using 6R robot type MH-12 YASKAWA, with a Table 3. Models of the 6R robot joints and the corresponding compliance matrices.
where V force is the 6 × 1 vector of three linear forces and three angular moments acting on the tool tip. V dev is the 6 × 1 vector of three linear and three rotational displacements of the tool tip.
The various components of the full deviation vector are calculated and analyzed. On the right side of the manipulation display screen in Figure 4   Let the WP plate be deflected at angle α from the vertical. XYZ is the coordinate system associated with the robot ( Figure 5); X'Y'Z' is the coordinate system associated with the plane of the WP plate. vDev = {δx, δy, δz} is the previously calculated vector of full deviation and its known components in the XYZ coordinate system, and vDev = {δx', δy', δz'} is the same vector in the X'Y'Z' coordinate system. The problem to find δx', δy', δz' and vector devPlane = {δy', δz'}.
[ ] The following four cases (a-d) of milling forces were considered for calculation of the deviations. Values of forces were estimated during measurements, described in Section 4.2.

Visualization the Workspace According to Limits of Deviations
The value of full deviation vector vDev as well as values of vectors of plane deviations for the devPlane was calculated for each FP. The value of the deviation at each FP was taken equal to the maximum value of the FP corresponding to all four cases mentioned above. The contour maps in Figure 6 identify regions where deviations do not exceed the required level. Figure 6(A) shows a complete contour map of plane deviations in the workspace of interest. Figure 6(B) shows the permissible part of the space where deviation limits do not exceed the predefined limits.
A set of contour maps for different distances between the robot's Z-axis and the central point of the working plane given in Table 4. Central point coordinates are {800, 0, 800}. The last third column, the permissible zones where deviations do not exceed the limit. The set of these figures can be interpreted as a set of cross-sections of the three-dimensional region of the permissible FP positions in terms of the allowable robot machining error.
Analogous set of contour maps for different angles between the working plane and the vertical (the robot's Z-axis) shown in Table 5.

Experimental Investigations with 6R Robot Type YASKAWA MH-12
The experimental part of the study consisted of milling a circular groove, on a flat WP surface, with 100-mm radius, using a high-speed mini spindle mounted on MH-12 YASKAWA robot. The considered sources of errors were the robot's positioning accuracy, which determines the specified trajectory of the tool motion, and TCP deviations from the specified path due to the finite rigidity of the robot joints. The obtained machining accuracy was compared with the calculated theoretical value. The effect of robot rigidity on the abovementioned deviations was determined. Milling process forces considered in the described study and within the developed model and it was therefore necessary to measure them.

Equipment
The milling experiment used the following equipment:  Robot MH-12 YASKAWA (1), spindle adaptor (2), (Figure 7(A)).  The high-speed Nakanishi milling spindle (3), (Figure 7(B)), which has a rotation speed of up to 60,000 rpm, a milling cutter with 2 mm in diameter.  The milled WP was fixed on the 3-component dynamometer as shown in Figure 7(A). The robot was programmed such that, the TCP moved along the same circular path as the ball-bar test path. Coordinates of the circle center are {1000, 0, 800} in robot coordinate system. A predicted deviation of less than 0.07 mm accrued during the abovementioned circular milling path ( Figure 8). The cutter feed along the path was 100 mm/min and spindle speed was 50,000 rpm. Ten milling cycles were performed with a depth of cut of 0.2 mm (total groove depth was 2 mm).

Description of the Experiment
The following steps were performed:  The robot was programmed to move the tool along a circular path.  Two full ball-bar tests were performed. The test results were combined into one chart and shown in Figure 11.  The forces acting on the tool were measured in three directions along the X, Y, and Z-axes simultaneously with the milling process. The measurement results were recorded and displayed in the form of graphs ( Figure 12) on the monitor of the Kistler measuring system.  After the milling path was completed, the precise values of the radius of the circle were measured every 5 degrees to define the groove roundness ( Figure  13).

Comparison and Analysis of the Experimental Results
To analyze the positioning accuracy of the robot without an external load, the measurements of the ball-bar tests were combined into one diagram, shown in Figure 11 and were compared.
The solid and dashed blue lines in Figure 11 represent the first and second ball-bar clockwise (CW) measurement series, while the solid and dashed green lines represent are the first and second ball-bar counterclockwise (CCW) measurement series. The difference between the measured values in CW versus CCW tests (not more than 100 μm) was much greater than the measurement differences found for the unidirectional tests: CW1 versus CW2 and CCW1 versus CCW2 (approximately 20 μm). Therefore, the main source of TCP position  The stiffness/compliance model of the MH-12 robot previously described in Section 2 was used to calculate the errors caused by the robot joint stiffness during circular motion. Figure 15 shows the pose of robot during milling the circular path. Engineering    Figure 16. The calculated position errors in the WP plane lie within a ±60 μm range for milling with the manifested conditions described in Section 4.2.
Precise distances from the obtained circular groove center to the planed groove centerline were measured, applying steps of 5˚ along the planed groove centerline as shown in Figure 13. As clear from Figure 13 the real deviations of the milled circle from the nominal one lie within the ±200 μm range. The calculated deviations due to stiffness of the robot joints were ±60 μm. The difference of ±140 μm was caused by the robot's positioning error.

Conclusions
The achievable accuracy and stiffness of serial-kinematic IRs used for machining tasks were investigated and experimentally verified. An actual machining process that is typically used for accuracy testing has been selected: profile milling of the circumferential slot by an end mill tool. High-speed machining with 50,000 rpm under optimal cutting conditions was performed.
To obtain maximally achievable IR machining accuracy, 1) a two-step teaching procedure was applied, it included point position correction using Renishaw ball-bar measurements and a final ball-bar test to establish the TCP motion trajectory; 2) a simulation tool based on the FSF approach, previously applied only for machine tools, to predict accuracy and stiffness within the IR workspace; and 3) the developed software tool, which allowed calculation and visualization of the compliance matrix for an arbitrary robot posture. Also, the given compliance matrix, joint rotation angles, and the values of the tool tip deviation vector components, caused due to the presence of cutting forces, were displayed. Applying the mentioned above procedure enabled us to predict and optimize the positioning error at a given point in the workspace.
The measured forces during the machining experiments were used to calculate deviations from the nominal position that occurred due to robot rigidity. The obtained deviation values on the selected machining plane were represented as contour maps. Different tilting angles of the WP plane were represented in a contour map. Such mapping makes it possible to optimize the WP position and orientation relative to the robot with respect to machining accuracy.
The obtained roundness of the milled circle was ±200 μm, which coincides well with the sum of deviations from the nominal radius that the robot positioning error (≈140 μm) and robot rigidity (≈60 μm) caused. Thus, the effect of positioning error on milling accuracy is 2.2 times greater than that of the robot joint stiffness. The proposed simulation tool is convenient for estimating robot positioning error due to robot compliance within a workspace with preliminary limited errors.