A Robotic Milling System Based on 3D Point Cloud

: Industrial robots have advantages in the processing of large-scale components in the aerospace industry. Compared to CNC machine tools, robot arms are cheaper and easier to deploy. However, due to the poor consistency of incoming materials, large-scale and lightweight components make it difﬁcult to automate robotic machining. In addition, the stiffness of the tandem structure is quite low. Therefore, the stability of the milling process is always a concern. In this paper, the robotic milling research is carried out for the welding pre-processing technology of large-scale components. In order to realize the automatic production of low-conformity parts, the on-site measurement–planning–processing method is adopted with the laser proﬁler. On the one hand, the laser proﬁler hand–eye calibration method is optimized to improve the measurement accuracy. On the other hand, the stiffness of the robot’s processing posture is optimized, combined with the angle of the ﬁxture turntable. Finally, the experiment shows the feasibility of the on-site measurement–planning–processing method and veriﬁes the correctness of the stiffness model.


Introduction
Large-scale components are commonly found in the aerospace and shipbuilding fields. The manufacturing requirements of these components require large-scale machine tools. However, large-scale or special-purpose machine tools are expensive. The deployment of large components on the machine tool is time-consuming and laborious [1,2]. Therefore, it is difficult to meet the demand. In addition, some workpieces have to use machine tools because manual work cannot guarantee product quality, though they do not have high accuracy requirements. Therefore, industrial robots are widely used in this case. Industrial robots have incomparable flexibility and price advantages to CNC machine tools [3][4][5]. They are easy to deploy and have a large working space. They have gradually become key pieces of equipment in the field of aerospace manufacturing. However, their shortcomings are also obvious. For example, the stiffness of the tandem structure is poor, the programming difficulty is often higher than that of the CNC system, and it is not easy to automate the incoming materials with low consistency. Especially in the milling process, the serial structure of industrial robots makes machine chatter unavoidable. Thus, the chatter problem has always been a research hotspot [6][7][8][9].
This paper focuses on the robotic processing of large structural parts in the aerospace field. Since large-size components are often welded by multiple parts, bevel processing is indispensable in the welding process. However, the bevel processing is often carried out manually. Because of its large amount of removal, workers need to finish this part by working in a dusty and noisy environment for hours. In addition, due to the low stiffness of this lightweight component, the consistency of its previous processing is very poor. Furthermore, the workpiece has no positioning reference, which brings great difficulties to automation. In this paper, the authors will accomplish robotic bevel processing for lowconsistency workpieces automatically. At the same time, considering that the processing method is milling, the stiffness performance of the robot needs to be optimized to improve processing stability.
Due to the structural characteristics and manufacturing process characteristics of large lightweight components, the incoming consistency of the bevel processing process is low. This, together with no benchmark fixtures, results in barriers to deploying automated production. Therefore, the robot cannot process with an inflexible program such as traditional CNC milling.
In this paper, a flexible robotic milling process framework is proposed, as shown in Figure 1. The robot carries a laser line sensor to scan and obtain the point cloud of the workpiece profile. As a laser profiler can provide high-accuracy data, the point cloud data can be processed and applied to trajectory planning directly. According to bevel processing technical requirements, the 3D cameras are not suitable for milling applications because of their low accuracy. However, the hand-eye calibration method of the profiler based on a laser line is not widely applied. To improve detection accuracy, a hand-eye calibration method for a laser line sensor based on Genetic Algorithms (GA) is proposed.
Compared to CNC, the disadvantages of robotic milling are obvious. The low stiffness aggravates chatter during the milling process. Thus, in this section, the effect of stiffness on processing is considered. Robotic milling stability is improved by optimizing the manipulator's processing posture. On the other hand, the velocity of the manipulator's trajectory is also optimized to improve the processing efficiency.
The remainder of the paper is organized as follows: The related works of several key techniques are introduced in Section 2. Section 3 will present a novel hand-eye calibration method. Section 4 will provide the feature analysis method and path planning based on point cloud data of the actual part. The milling postures will be optimized in Section 5 based on stiffness modeling. The results and discussion of the system experiment will be shown in Section 6. Section 7 will conclude the paper.

Related Works
This section includes the robot hand-eye calibration method for a laser profiler and research on the stiffness of robots.
Aerospace structural parts often require online measurements to obtain the structural characteristics of the current workpiece. According to different applications, the accuracy requirements are different, and the sensors are also different. In the existing research, there are many methods to measure the whole workpiece or the region of interest using robots combined with measurement sensors. Kuss et al. used a 3D camera to achieve a three-dimensional measurement of structural steel workpieces [10] and completed the deburring process according to the measurement results. Han et al. used a system of robots and 3D cameras to measure the three-dimensional shape of the workpiece on the forging production line [11]. Yu designed the line laser rotating mirror scanning equipment and successfully applied it to the scanning of the car body in white by the robot system [12].
Wang et al. used a mobile robot combined with a 3D camera to increase the scanning range of the robot and realized the scanning of the wind turbine blade model [13]. Ge et al. used robots and line laser sensors to realize the welding seam of the steel pipe workpiece via online measurement and removal [14]. Guo et al. used robots and line laser sensors to measure and remove spiral welds on steel pipes [15]. Among them, it can be found that due to its limited accuracy, 3D cameras are often used in applications where precision requirements such as deburring and polishing are not high. Sensors based on the principle of laser imaging, usually used for quality inspection or welding seam processing, requires high precision or simple scanning methods. In this paper, for bevel processing, the accuracy requirements are relatively high, and 3D cameras are not suitable. Therefore, the laser profile sensor is the first choice.
Generally speaking, in order to adapt to the task of 3D measurement of large-scale components, the system adopts the form of "eye on hand". Hand-eye calibration methods currently mainly include traditional methods, optimization methods, and some probabilitybased methods, but most of them are for camera sensors. The research ideas of hand-eye calibration for a laser profiler are similar and simple. A typical calibration method is to use a robot to measure a standard ball in multiple poses [16][17][18][19] and construct a similar equation to solve it. In addition, some scholars use some special-shaped calibration objects, such as disc-shaped [20], X-shaped [21] and multi-step calibration objects [22]. Carlson et al. used a three-orthogonal plane method for calibration [23], and Sharifzadeh et al. improved this by using a single plane to achieve convenient and stable calibration [24]. However, for the measurement task, what needs to be obtained is the 3D reconstruction result of the measured object. From the perspective of 3D reconstruction, there are many systematic errors, including distortion and scattering, which cannot be solved by traditional methods. Therefore, on this basis, optimization methods may be more suitable to improve data accuracy.
Industrial robots are composed of tandem rotary joints. Although this structure is beneficial to increase the working space of the robot, it will result in significantly low stiffness characteristics. The cutting force acting on the end of the robot arm during the machining process will cause undesired displacements. Due to the low stiffness, the error produced during the machining process accounts for the main proportion of the machining error [25]. Therefore, a lot of research has focused on the stiffness modeling and stiffness error compensation of the robot. Stiffness modeling methods can be divided into two categories: virtual joint method (VJM), which describes elastic components as lumped parameter models [26,27], and finite element analysis (FEA), which calculates elastic deformation based on computer finite element-aided design tools [28][29][30]. The virtual joint method can easily obtain the stiffness model of the robot, but the accuracy is not as good as the finite element method. The main advantage of the finite element analysis method is its higher accuracy, but it requires higher computing power support and higher modeling skills. Gong [31] believes that the influence of robot load and the weight of the connecting rod on the deformation of the joint is significantly greater than the flexible deformation of the connecting rod, so most of the literature only considers the influence of the flexible joint deformation [31][32][33]. Khalil [34] uses the cantilever beam model to consider the flexible deformation of the connecting rod, which is too complex to be suitable for engineering applications. Therefore, the virtual joint method is widely used in robot processing applications because of its high computational efficiency and acceptable accuracy. Salisbury [35] first proposed a robot stiffness model. The transmission components of each joint of the robot (such as reducer, actuator, etc.) are the main sources of flexible deformation of the robot. A virtual torsion spring is used to represent the stiffness of each joint, thereby establishing a classic stiffness model. However, this model is only valid when the robot has no load. Therefore, Chen et al. proposed a Conservative Congruential Transformation (CCT) [36,37] . CCT describes the relationship between the Cartesian stiffness matrix and the joint stiffness matrix and considers the displacement change of the robot caused by the load on the end of the robot. This paper is based on this method to identify the stiffness of the robot. The robot processing path is optimized using this stiffness index.

Modeling
The robot hand-eye calibration is used to obtain the coordinate system of the sensor relative to the robot end. The calibration model of the robot line laser scanning system is shown in Figure 2. In order to calibrate the hand-eye relationship, the following coordinate system needs to be established: {C 0 } : o 0 − x 0 y 0 z 0 stands for robot base coordinate system; {C e } : o e − x e y e z e stands for robot tool coordinate system, where z e is the axe of robot flange center; {C c } : o c − x c y c z c stands for the profiler coordinate system, where y c is perpendicular to the laser plane. It is easy to obtain the transformation relationship T e 0 through the kinematics between {C e } and {C 0 }. Therefore, the goal is to solve the transformation relationship T c e between {C c } and {C e }. Let P c ball be the center of calibration ball in {C 0 }, and P 0 ball be that in {C 0 }. Therefore, the relationship is: Equation (1) can be transformed into the following form: where r and t denote the rotation and translation of T. The relative relationship between the center of the laser calibration sphere and the profiler coordinate system is shown in Figure 3. Therefore, according to the measurement principle of the line laser sensor, we use the least square method to solve the radius R plane and the secant plane center O n (only x and z of P c ball ) of the light plane and the laser calibration sphere. Thus, h (y of P c ball ) in Figure 4 can be calculated with R ball and R plane .  As the calibration sphere is fixed, P 0 ball is constant in {C 0 }. To solve T c e , equations can be set from different positions and postures of the manipulator, that is: · · · P 0 ball = n r e 0 (r c e × n P c ball + t c e ) + n t e 0 (3)

Solution
As the analytical method is applied, the results show a great degree of dispersion in y values of P c ball . Therefore, a solution based on a Genetic Algorithm (GA) is proposed for the calibration method to improve the accuracy. GA is an optimization process that is suitable for a wide range of values, many parameters and nonlinearity. Since directly solving the six parameters of the hand-eye relationship is easy to converge to the local optimum, this part adopts the method of applying GA to solve the rotation component and translation component, respectively.
As the initial value of the rotation and translation components are calculated from the designed model, the objective function and individual fitness functions are designed as follows. The rotation matrix of the x component is: Then the rotation of T c e is: Thus, the objective function can be expressed as: where · 2 denotes the Euclidean norm, and roll, pitch, yaw denote Euler angle inputs. Floating point encoding is adopted instead of binary encoding to shorten the encoding length and improve the efficiency. The individual fitness function can be expressed as: Let t c e = [a, b, c] T , while r c e is calculated above. The same is true for the translation component; the objective function can be expressed as:

Calibration Results
To obtain a reasonable calibration results, 3 time measurements are carried out for each of the 3 postures. The analytical method and GA method are both tried out, so there are 18 groups of data, as shown in Table 1. In Figure 5, calculated sphere centers are shown in an intuitive form. The standard deviation of the sphere centers by analytical method is σ x = 0.6793, σ y = 3.3567, σ z = 0.8215, while that by the GA method is σ x = 0.0660, σ y = 0.4136, σ z = 0.1131. It shows that the method based on GA has a good robustness to system errors. The final hand-eye calibration result is t x = 126.92 mm, t y = 9.44 mm,  The reason for the inaccuracy of the analytical method is that the diameter of the calibration sphere used in this calibration is a bit small (10 mm), while the repeated positioning accuracy of the robot is 0.2 mm. The sensor's sensitivity to the light environment and the accumulation of errors, etc., lead to a large deviation between the results of this calibration using the analytical method and the design value. According to the above factors, the authors design the corresponding objective function for the GA on the basis of the analytical method and obtain a more accurate result.

Scanning
The workpiece profile in point cloud data is the precondition of all the analysis jobs. The manipulator carried a profiler to scan the workpiece, which is shown in Figure 6. The main problem is how to align the data obtained from different devices. Specifically, two kinds of data are needed to reconstruct the point cloud data, tool position information from the manipulator and point information from the sensor. For every profile points group P ic , there is a tool position T e 0 . Therefore, multiple measurements can be expressed as: where T c e denotes the hand-eye calibration matrix in the previous section, and n denotes the times of measurements. Then, the final point cloud data are fused as However, to reconstruct a high accuracy workpiece model, n needs to be large. Multiple measurements will consume a lot of time. According to the profiler's communication performance, the data frequency can be 50Hz. Continuous measurement by a moving manipulator can solve this problem. The manipulator is configured to move in a straight line (easy to interpolate) at a uniform speed. The authors use time stamps to align the sensor data with 1 T e 0 and n T e 0 . Thus, other T e 0 can be calculated by interpolation. Then the final point cloud data are computed with Equation (9). Therefore, this method can improve the efficiency significantly. Figure 7 shows the situation of robotic milling on a workpiece. As the process features are suitable for peripheral milling, the end side face of the workpiece needs to be extracted for tool path planning. In the previous part, the scanned model is made up of 2D profiles. Every profile is one frame, and the edge points are made of vertexes in every frame. As vertex is extremum, it is easy to extract in a 2D plane. Figure 8a,b shows the performance of edge extraction.  The least square method, a simple identification method, is used for plane fitting. From the plane expression AX + BY + C = Z, the expression in matrix AX = b form is:

Extraction
The objective function is: Figure 8c shows the results of plane extraction. Preparing for tool planning, the plane needs dimension reduction, which can simplify the tool path generation, as shown in

Planning
When considering the cutting amount, the z-axis of the peripheral milling cutter should be planned on the direction of a narrow feather. Then the cutter motion direction y-axis is perpendicular to the z-axis and normal to the extraction plane x-axis. After the coordinate is defined, a series of key points of the tool path can be generated every 1 mm, as shown in Figure 8e.
Finally, to face the requirements of processing, the angle of the edge side face and bevel plane is 40 • , and the tool path also needs to be rotated around the z-axis by 40 • . Thus one tool path is generated, as shown in Figure 8f, and all other tool paths are spatially translated along the edge side face.

Processing Posture Optimization
Mode coupling chatter was identified as the dominant source of vibrations in robotic machining, largely due to the inherent low structure stiffness of an industrial robot. In this section, the authors establish the stiffness model and identify the stiffness matrix. Then, a new evaluation indicator is proposed and applied to optimize the manipulator's stiffness performance.
Before modeling, three hypotheses should be stated here: (1) The end effector is infinitely stiff. (2) Links of the manipulator are infinitely stiff. (3) All the processing is in a steady state, as the feeding motion is quite slow. Therefore, all the deformations are caused by joint flexibility.

Cartesian Stiffness Model
In the process of robot grasping objects and operations, K x is regarded as a measure between the displacement of the robot tool coordinate system and the interaction force when an external force acts on the robot end effector. In the classic stiffness description, it is related to the joint stiffness of robot K θ as follows: According to the Taylor expansion formula, the relationship between the small force vector d f and small displacement vector dx on the end of the manipulator is expressed as (higher-order terms are ignored): According to the robot Jacobian matrix, dx can be expressed as: However, the classic stiffness description is valid only for a quasi-static robot configuration with no loading. During processing, the milling force can cause joint torque. As joint torque dT can be given by J T K x dx, Equation (13) is substituted into a virtual work principle, which is written as: The last term can be simplified by using the method proposed by Chen et al., that is: By combining Equations (14)-(16), Cartesian stiffness K x is given as: Thus, the additional stiffness term due to external loading has been incorporated using the Conservative Congruence Transformation (CCT) approach.

Joint Stiffness Identification
Judging from Equation (17), the joint stiffness is very sensitive to the condition number of the Jacobian matrix. In order to ensure that the stiffness identification program can converge to a stable value, it is necessary to determine the dexterity in each posture to select a posture suitable for stiffness identification.
The physical meaning of the dexterity of a robot is a comprehensive measure of the movement ability in all directions, which is used to measure the flexibility of the robot in various spatial poses. It is defined as: where λ i are the eigenvalues of J J T , and σ i are the singular values of the Jacobian matrix.
As the dexterity is deeply affected by Axis II and Axis III of the 6 DOF robot, a contour map of ω −1 under the two axes is given in Figure 9 to analyze the suitable postures. In addition, Axis IV, V, VI are configured at 45 • to avoid singular points. To simplify the identification, the authors hope to ignore K c term. Equation (17) shows that the higher the milling force is, the bigger value the K c gets. A 2000 N force vector and a 200 Nm torque are applied to the end of the manipulator. Here, two indices f p and f r are adopted to describe the influence of K c on K x based on translational and rotational displacements. They are defined as follows: and where δ p Kc and δ pK c denotes the point-displacements of robot end-effector under milling force with and without considering K c , respectively. In addition, δ r xKc , δ r yKc , δ r zKc and δ rx Kc , δ rz Kc , δ rz Kc denote the small rotations of the robot end-effector under torque with and without considering K c . Figure 10 shows the contour map of f p and f r on Axis II and Axis III. The darker color means the lower influence. However, all the values of f p and f r are tiny. When f p < 0.01 and f r < 0.0001 rad, the K c term is considered negligible. Therefore, substitute Equation (12) into Equation (13), the force can be written as: or If joint flexibility is expressed as: Then, by expanding Equation (22), it turns out that: Consequently, Equation (24) is the form for force loading tests. According to the analysis above, a posture of higher dexterity and lower influence of K c is chosen. Then, the joint stiffness can be computed by the amount of applied force values and displacements of the robot's end-effector. In order to obtain a higher accuracy stiffness matrix of joints, the least squares method will be adopted with multiple measurements. The identification experiment is carried out in Section 6.2.

Optimization
Suppose the stiffness matrix of joints has been obtained, Cartesian stiffness can be calculated at every posture by Equation (14), thus providing a different measure for the influence of milling force at different postures. Figure 11 shows the robotic milling system, including an industrial robot, an electric spindle (end-effector), and a fixture turntable. In this system, there are two redundant degrees of freedom, rotation of the fixture turntable and rotation of the cutter. Therefore, for the milling job, there are always a series of processing postures gaining better stiffness performance. In addition, two conditions must be met: (1) the postures of the manipulator are solvable, (2) the postures are non-collision.
An optimization problem is defined as follows: where θ i denotes the rotation offset of the milling cutter and α i denotes the rotation of the fixture turntable, i denotes the points in the whole tool paths, and Θ denotes 6 angles of the posture in joint space. The objective function is the average of the displacements δx, which are calculated according to the stiffness equation. Figure 11. The robotic milling system.

Experiments, Results and Discussions
In this section, simulations and experiments are carried out to verify the methods proposed above. This part contains the complete process of the robotic milling system, including measurement and planning, stiffness identification and optimization, and machining test.

Tool Path Planning
The coordinate system is defined as shown in Figure 12. O tool is the origin of the robot tool coordinate system. O p is the origin of the processing coordinate system. The industrial robot is ABB IRB6700-150 with a 150kg payload. Based on the single up-milling tool path generated in Section 4, all the tool paths are planned layer by layer, as shown in Figure 13. However, not all tool paths cut the workpiece. Considering the processing efficiency, the authors accelerate the motion of the non-cutting path. A velocity planning method is proposed here: where d min denotes the minimum distance between the y-axis of the processing coordinate system and the point cloud of the workpiece; v f ast and v slow denote the velocity of the cutter at non-cutting points; R denotes the radius of the cutter. Figure 14 shows the performance of the non-cutting path velocity planning. The darker color denotes a lower velocity and lighter denotes higher.   Figure 15 shows the platform of the stiffness identification. The authors use a tracker to measure the displacements of the end-effector. A tool is designed for the installation of the target ball and force sensor. A bench vise is adopted to load the force. According to Section 5, three indices need to be considered, namely the dexterity, f p and f r , shown in Figure 16. Three points, denoting three postures, are chosen from the three charts to carry out the identification experiment. The detailed data are partly shown in the Appendix A. The 6D displacement vector can only be solved accurately when the end of the robot is in micro-motion. Only in this case, the relationship formula between the stiffness of the robot joint and the end of the Cartesian stiffness is valid. Therefore, the maximum load during the experiment is set to 200 N.

Stiffness Identification and Optimization
As the installation location of the sensor does not coincide with tool coordinate system, transformations of the force and torque are necessary. For force, the transformation is: where f f is the force data of the sensor, and e f is the force of the end of the robot. e T f denotes the transformation matrix from the robot to the sensor. For torque, that is: where τ denotes the torque on the end of the robot. δ denotes the distance from the sensor to the end of the robot. e R tool denotes the rotation matrix from the end of the robot to the tool coordinate system. Finally, the least squares method is applied to multiple pairs of force and displacement to improve the accuracy. The detailed data of robot posture configuration and loading are presented in the Appendix A. The joint stiffness of the robot is shown in the following table (Table 2). To optimize the angles of the two redundant DOFs of the machining system, the function in (25) is defined as the stiffness evaluation index. Figure 17 shows the result of the traversal simulation of the two redundant DOFs of the machining system using the stiffness evaluation index. Figure 17b is the three-dimensional figure of the simulation result, and Figure 17a is the contour sketch of the right figure. Figure 17c,d is the robot postures before and after optimization. The technological parameters of the milling process: the spindle speed is 6000 r/min, the feed speed is 1 mm/s, the feed is 1mm, and the milling cutter is a 4-blade end mill with a diameter of 10mm. The milling force simulation program and the method of milling force parameter identification refer to the method introduced in [33] for simulation calculation.

Machining Test
It should be noted that the maximum mean square error of the hand-eye calibration reprojection is 0.41 mm. According to the analysis and optimization results of the processing accuracy in the previous part, the absolute positioning accuracy of the robot improves to be ±0.1 mm by using the hand-eye calibration method. The optimized processing gains the average displacement mode length of the process at 0.07 mm. Considering the above data comprehensively, it is not difficult to find that the machining results cannot meet most of the needs by relying on the absolute positioning accuracy of the robot for milling processing. Therefore, manual tool setting is required before milling to minimize the impact of the absolute positioning accuracy of the robot on the machining accuracy.
For the robot bevel milling process, the comparison effect of the processing result and manual bevel processing are shown in Figure 18. It can be found that the surface quality under visual inspection is similar between robotic milling and manual grinding. The bevel inclination angle needs to be obtained by applying numerical software analysis to the point cloud data after the bevel milling is completed, and the average distance of the blunt edge also needs to be obtained by numerical analysis software. Before evaluating the processing result, the processing requirements are shown in Figure 19.  The distance between the upper and lower curves on the green plane is the required root size. The result of the distance between the upper and lower edge of the bevel end surface obtained by point cloud data measurement is shown in Figure 20. The basic size is 1.5 mm. The shadow is the tolerance zone in the process requirements. It can be seen that after removing the influence of some error points, the overall indicators are within the process requirements.  Figure 21 is a schematic diagram of calculating the inclination angle of the bevel based on the workpiece point cloud after processing. Nine reference planes are selected equidistantly on the bottom edge curve. Based on numerical analysis software, the bevel angle in each section is analyzed and marked. The error of the angle is under ±1 • , which meets the processing requirements.

Conclusions
In this paper, the authors present a completed work on a robotic milling system. The significance of the work is that it can realize automated robotic milling of lowconsistency parts. It liberates workers from the harsh environment and improves the production efficiency.
The main contributions of this paper are as follows: (1) A hand-eye calibration optimization method is proposed based on GA, which presented a good robustness to system errors.
(2) The authors propose a processing feature extraction method based on the point cloud data of the actual workpiece. The tool paths are generated based on velocity optimization.
(3) A processing posture optimization model based on the robot stiffness is proposed for the specific system with a fixture turntable. The stiffness model ignoring the Kc term is verified by simulation.
(4) The robotic bevel milling task is completed, and the results meet the processing requirements.
In future work, the authors will conduct more process tests to optimize process parameters. The milling path planning method proposed in this paper has good regional adaptability to the different postures of the workpiece, but the robotic milling is not limited to the workpiece involved in this paper. The follow-up should conduct detailed research on the regional features of other parts to improve the adaptability of the tool path planning method.

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