Novel surface optimization for trajectory reconstruction in industrial robot tasks

This study presents a novel algorithm implementation that optimizes manually recorded toolpaths with the use of a 3D-workpiece model to reduce manual error induced. The novel algorithm has three steps: workpiece declaration, manual toolpath declaration, and toolpath optimization using steepest descent algorithm. Steepest descent finds the surface route wherein the manually recorded toolpaths traverse over a 3D-workpiece surface. The optimized toolpaths were simulated and tested with an industrial robot showing minimal error compared to the desired optimized toolpaths. The results obtained from the presented implementation on three different trajectories demonstrate that the proposed methodology can reduce the manual error induced using as a reference the CAD-workpiece surface.


Introduction
The industrial robot-workpiece interaction has an important role to perform tasks over complex surfaces such as painting 1-3 milling, 4 grinding, 5 laser cladding, 6 nondestructive examinations, 7,8 and robotic automated fiber placement applications. 9,10 The path planning can be classified into online and offline programming. In online programming, the main disadvantages are the time spent on manual robot programming and the manual error induced (error due to operator skills-dependency). On another side, the performance of off-line path planning methods is faster and often model-based 3D CAD. 7,[11][12][13] Some robotic applications take advantage of 3D-workpiece for diverse applications, for example, in Xingguo and Li, 14 a robotic model-based offline programming method made the projection of a 2D-plane onto the surface of a CAD model for machining processes, taking into consideration the spin angle of the tool. If a CAD model is not available, it is possible to take advantage of reverse engineering techniques using scanned data 15,16 offering discrete 3Dworkpiece models where special algorithms generate trajectories on point clouds and meshes. 8,13 In an earlier study, 3 the authors propose an incremental approach for automated trajectory generation from CAD model, the whole surface is meshed and divided into surface sections (SS) based on spray radius, and these SS are painted by spray gun located at a standoff distance from the surface section point (SSP) along the surface section normal (SSN). To calculate the SSP, the weighted average of all the centroids of the triangle inside the SS is used and the weighted average of the normals of the triangles in the SS determines the SSN. Subsequently, a new SS is generated with an incremental distance along the x axis, steps are repeated till the end of the surface is reached, and the incremental distance in y axis depends on the optimal overlap distance. To generate the optimal overlap distance and the velocity for a given paint pass (SS), a genetic algorithm is used and the process ends when the whole surface is reached. Li et al. 6 propose a path planning method whose key lies in the topological reconstruction of the artificial joint surface. From an STL file, the information to create the relations among points, edges, and facets is extracted, facilitating the path planning method. Via equidistant parallel planes that intersect the CAD model, a series of intersections point were generated (slicing algorithm) to obtain transversal lines that use the topological relation of the edges instead of using a sorting algorithm to reduce the number of intersection calculations. The height error method is used for the robot interpolation points and the normal vectors of each point can calculate the position and pose of the robot tool center point (TCP) to ensure the perpendicularity between the laser beam and the surface. Finally, by connecting the interpolation points sequentially, the motion path is obtained.
Other works have addressed the problem by mesh following techniques using a triangular mesh of a surface as a guide for path-planning algorithms; Mineo et al. 7,8 present a flexible off-line trajectory planning for inspection of complex curved surfaces on nondestructive testing systems. The toolpath is generated to follow the contour of a triangular meshed CAD surface without the need for an approximating analytical surface. To generate the path planning, it is necessary to create a curve contour of the surface, as the surface is a triangular mesh, the surface edge is a segment with extremities of two adjacent triangles, then always is possible to find the intersections points between the plane perpendicular to the segment and the edges of the triangles of the surface mesh. The curvilinear distance is calculated by cumulating the segment of the intersection points from the reference edge that is the distance along the surface contour where every point of the curve must be normal to the relative triangle.
This article presents a novel algorithm implementation that combines the skills of the operator and the advantages of the CAD model to optimize manually recorded toolpaths with the use of 3D-workpiece models to reduce manual error induced avoiding exhaustive touch up reprogramming when considering large, workpiece surfaces. The procedure was tested both with a simulation and experimentation on a real industrial robot arm; model ABB IRB 1600-7/1.45 type A. The results show a high-quality trajectory reproduction, allowing more accurate adjustments on manually recorded toolpaths.
The organization of this article is divided into four sections: Methodology presents a generalized explanation of the proposed algorithm for optimizing manually recorded trajectories, Steepest Descent algorithm selects the number of nodes in the objective function for reconstructing the surface route on the 3D-workpiece surface. Experimentation and Results were carried out at ABB industrial robot using three different trajectories where the mean absolute error was calculated by comparing the robot TCP trajectory and the desired optimized trajectory. Finally, the Conclusions and Future Work are presented.

Methodology
The presented novel algorithm implementation combines the manually recorded toolpaths stored in a text file with the use of a 3D-workpiece-model to optimize robotic toolpaths and reduce manual error induced. The steepest descent algorithm finds the surface route wherein the manually recorded toolpaths traverse over the CAD-workpiece. The three-step procedure of the novel algorithm implementation is presented as follows: Step 1: Workpiece declaration The first step is the workpiece declaration, Figure 1 shows a curved surface that represents a 3D-workpiece using point cloud data. This figure includes 656 points, covering an area of 0:24 m 2 . The points were uniformly distributed over the surface but can also be randomly distributed. The application of the presented methodology focuses directly on surjective CAD surfaces on 3D as in previous works, 8,13 where every z coordinate in the codomain Z from the workpiece surface has the corresponding xy-couple coordinate such that z ¼ f ðx; yÞ. The shape of the CAD-workpiece surface and its coordinate reference frame determine if the workpiece is surjective. The manually recorded toolpaths and the point cloud data that conform to the CAD-workpiece can be referenced to the robot base coordinate system or any other specific coordinate system (see Figure 2). This characteristic allows identifying the distancing hx w , hy w , and hz w from any coordinate system with respect to any point from the CAD surface, wherein any point can be defined as (DBW i ) for i ¼ number of points that conforms the 3D-workpiece. The presented methodology also considers a TCP in the flange of the robot or a tool attached to it.

Step 2: Manual toolpath declaration
The second step is the declaration of manually recorded toolpaths. These manually recorded toolpaths conventionally present variable standoff distancing between TCP and CAD-workpiece (represented as point cloud data), due to manually induced error. Every manually recorded toolpath can be traversed using at least two poses sequentially, and the result can be a straight toolpath that is subdivided into tiny segments. 17 As an example, Figure 3 shows five poses ðP 1 ; P 2 ; P 3 ; P 4 ; P 5 Þ to create four straight toolpaths. 17 Every pose can be represented as shown in equation (1), where hx i f ; hy i f , and hz i f correspond to the TCP of the robot, and E i z ; E i y , and E i x represent the orientation in Euler angles, such that Euler ZY X ða i ; b i ; g i Þ, where a i ; b i ; and g i are the rotation about Z, Y, and X axes, respectively. 18 Two successive poses defined as P i gpp and P ðiþ1Þ gpp can be expressed as in equations (1) and (2). For an "n" number sequenced toolpaths, it is possible to define P gpp as a matrix of M nx6 and  Every element of P gpp can also be expressed by a matrix of 4 Â 4 elements, where the first three columns represent the orientation and the last column is the corresponding position; for this case, P i gpp is expressed as in equation (4), where sðdÞ ¼ sinðdÞ and cðdÞ ¼ cosðdÞ, as a simplification for any angle d As a result of implementing equation (4), 17 straight toolpaths subdivided into tiny segments were created as a first approximant reference of the desired optimized trajectories.

Step 3: Toolpath optimization
In the third step, every segment from the straight toolpaths "DP i;iþ1 gpp " is optimized using the steepest descent algorithm. The algorithm finds an approximation of the 3Dworkpiece surface, where every DP i;iþ1 gpp traverse over the CAD-workpiece. The vector DP i;iþ1 gpp represents every possible transition from P i gpp to P iþ1 gpp ; therefore, this term can be converted from Euler angles to a representation using homogeneous transformation matrices, DP i;iþ1 gpp , P i gpp , and P iþ1 gpp E M 4x4 are represented similarly to equation (4) 17 (see equation (5)). These matrices correspond to the translation and orientation with respect to the TCP. Thus, it can be concluded that any transition from P i gpp to P iþ1 gpp can be represented as in equation (5), where Dhz i;iþ1 f corresponds to the initial estimation x ð0Þ , of the design variable x for the optimization problem (see Figure 4) The first two elements Dhx i;iþ1 f and Dhy i;iþ1 f in equation (5) are kept constant, whereas the third element Dhz i;iþ1 f forms part of the objective function as the initial estimation x ð0Þ of the design variable x. This means that objective function converges to an optimal by moving along the Z axis from the initial estimation Dhz i;iþ1 f as x ð0Þ with respect to the 3D-workpiece surface. 13 For every DP i;iþ1 gpp and Dhz i;iþ1 f as x ð0Þ that traverse over the CAD-workpiece (see Figure 5(a)), it is necessary to find a "j" number of points defined as DBW j to calculate the desired optimized surface route. These points DBW j can be found calculating the Euclidean norm as the difference of f as x ð0Þ . A priority queue is ordered using the Bubble sort algorithm and a radius of neighborhood "r À neighborhood," where "j" contained in DBW j is any number from "w" points chosen to optimize each segment of DP i;iþ1 gpp . Since xy coordinates are not modified for every DP i;iþ1 gpp , only DBW j 3;1 (equation (6)) is evaluated in the object function The distance between every transition d ss is defined as the difference in between DP i;iþ1 gpp and DP iþ1;iþ2 gpp and depends on the task and the desired accuracy, very small point intervals lead to a time increase in machining processes, whereas longer distances generate a rougher surface 19 (see Figure 5(b)). For continuous sequences of manually recorded poses such as P i gpp , P iþ1 gpp , P iþ2 gpp , . . . the linear path-planning correlates d ss with the sum of all the distances that conform a continuous path, P AB (see equation (7)), where t inc indicates normalized increments from 0 to 1, the normalized time interval Y represents all possible segment transitions on P i gpp , P iþ1 gpp , P iþ2 gpp , . . . 13 Steepest descent algorithm implementation The steepest descent (SD) algorithm for unconstrained optimization problems is a popular method for minimizing differentiable functions. 20 For this case, every segment DP i;iþ1 gpp is individually analyzed using DBW j 3;1 and SD algorithm. While Dhx i;iþ1 f and Dhy i;iþ1 f from every DP i;iþ1 gpp remains constant, Dhz i;iþ1 f moves along to Z-axis to approximate to the surface (see Figure 6). Equation (8) presents the objective function that minimizes the sum of the squared residuals (where a residual is defined as the difference between the chosen points from the 3D workpiece surface DBW j 3;1 and the design variable x). The number of terms in the objective function w increases as more elements of DBW j 3;1 are selected to reconstruct the surface. Since this is a quadratic function, it is possible to calculate the step size with equation (9) Thus,ŝ ðkÞ is obtained from f ðxÞ for the case where l k is dynamic as depicted in equation (10), while the approximation of x ðkþ1Þ is as in equation (11) x ðkþ1Þ ¼ x ðkÞ þ l kŝðkÞ A variant of steepest descent algorithm applied to the optimization of manually recorded toolpaths is presented in Algorithm 1. The stopping criterion for algorithm 1 is defined as df x ðkÞ =dx < e, where e is a constant that compares the change rate with respect to the objective function. Once every segment DP i;iþ1 gpp from every straight toolpath that has been optimized, the desired optimized toolpaths can be smoothed to avoid any possible discontinuity 21 (see Figure 7(a) and (b)). Finally, the desired optimized toolpaths can be displaced along the Z axis of its corresponding coordinate frame, with a standoff distance between the 3D-workpiece surface and the desired TCP (see Figure 7(c)).
Algorithm 2 shows the complete procedure used to optimize trajectories in the industrial robot for both the simulation and experimentation implemented in the following section.
find the Euclidian Norm l of DBW O and DP i;iþ1 gpp : execute the smoothing process, using DB 2 as input, based on a penalized least squares method, DCT, and IDCT; 25 execute the program corresponding to the robot kinematics (simulation on Matlab ® ); 26 execute the post-processor (for ABB ® systems) /* conversion of Matlab ® code into RAPID for experimentation 27 end

Experimentation
The experimentation was carried out using an ABB industrial robot IRB 1600-7/1.45 type A, (Figure 8) and a program-like based on algorithm 2 that converts a set of instructions in P gpp format into RAPID language (the language to program ABB robots). The performance of algorithm 2 was evaluated using three baseline trajectories over the 3D-workpiece surface. The first trajectory represents a closed trajectory with sharp corners, the second a zig-zag trajectory with sharp corners, and the third a closed trajectory with rounded corners [13]. Robot TCP motion was recorded and compared with desired trajectories to measure mean absolute error (MAE).
Trajectory 1 presents variable standoff distance in the Z axis simulating manual induced error during point-to-point programming or lead-through method (see Figure 9(a)). In this example, the trajectory finalizes its trajectory in the same pose as the beginning and it is composed of five poses as shown in Figures 3 and 9(a).
Trajectory 2 depicts a zigzag programmed in 2D (with components only in xy coordinates), where the trajectory has repetitive sharp corners (see Figure 9(b)). In this example, 14 poses were declared to create trajectory 2, wherein the start and end form an open trajectory.
Trajectory 3 presents a closed trajectory without sharp corners. The main purpose of this trajectory is to move the  TCP as close as possible from one corner to another without stopping. For this case, 41 poses characterize the trajectory, as presented in Figure 9(c).
The 3D-workpiece surface from Trajectories 1-3 (see Figure 9(a) to (c)) is known and was selected to validate algorithm 2 using an industrial robot. Equation (12) represents the 3D-workpiece surface with defined intervals of 0:6 hx w 0:9 and À0:4 hy w 0:4; hz w ¼ f ðhx w ; hy w Þ, where hx w ; hy w ; hz w can be any coordinate vector in the Euclidean space. The surface was represented as a point cloud, with three different densities 656, 2511, and 9821 points, corresponding to separations of ¼ 0:02; 0:01; and 0:005 m, no additional information from the 3D-workpiece was required f ðhx w ; hy w Þ ¼ 0:5 þ 1:8 Algorithm 2 was computed in MATLAB ® with a graphic user interface that includes the 3D-workpiece surface, robot, toolpath, and an additional window that allows selecting a file with P gpp poses (see Figure 10). The application simulates the robot motion performing the optimized desired trajectory before its practical implementation. At the end of the simulation, the application translates automatically the P gpp poses into RAPID code (the native language from ABB robots). The resulted file was loaded into the robot controller using a USB drive. The rapid code   created in MATLAB also includes additional instructions to store the real TCP location from the robot, via serial communication in real time using a cable that connects the controller with a PC.

Results
The three presented trajectories in the previous chapter were simulated and evaluated using the ABB Industrial robot IRB 1600-7/1.45 type A (Figure 8). Other configuration parameters are w for w ¼ 3 and the distancing defined as d ss for d ss ¼ 0:0025 m. Table 1 includes the number of poses for each trajectory before and after the implementation. This information serves to calculate the robot programming lines required. The number of segments generated in Table 1 depends mainly on d ss .
The mean absolute error (MAE) was measured from the difference between the desired trajectory route based on equation (12) and the TCP final measurement from the robot (reading the angular motion from encoders/resolvers and the length of every link). Table 2 shows the MAE obtained from the robot varying the mesh density (), for ¼ 0:02m; 0; :01m, and 0; :005m. Table 2 shows the MAE when ¼ 0:005 m, taking into consideration the individual average error for every reference coordinate in a trajectory and the total error as the sum of them.
The results obtained from Table 2 demonstrate that the mesh density () impacts directly on the reduction of the average error produced for the three trajectories from an average error of 0.4 mm to 0.1 mm. Likewise, Table 3 shows the MAE obtained from the same trajectories sorting the MAE by coordinates. The methodology presents comparable results with other methodologies taking into consideration the proposed mesh density 7,13 The mesh density and distance between every segment transition d ss directly increases or decreases the accuracy of the optimized trajectory. The error produced in xy coordinates was at some point expected and might be attributed to external factors from the proposed methodology such as robot serial configuration, low structural rigidity, and gear backlash 22

Conclusions and future work
In this article, a novel algorithm implementation combines the manually recorded toolpaths stored in a text file with the use of a CAD-workpiece model to reduce manual error induced and optimize the final robotic toolpath. The steepest descent algorithm finds the surface route wherein the manually recorded toolpaths traverse over the CADworkpiece.
The results shown in the previous chapter demonstrates that the proposed methodology can reduce the manual error induced using the CAD-workpiece as a reference. The results obtained for the three different trajectories evaluated indicate that the accuracy of the trajectory depends directly on the mesh-distancing , the distance between every segment transition d ss ¼ 0:0025 m, and external factors such as robot serial configuration, low structural rigidity, and gear backlash. The quality of the 3D-workpiece surface reconstruction also depends on the objective function. As for now, the implemented methodology only utilized three points from the 3D-workpiece surface for every DP i;iþ1 gpp and Dhz i;iþ1 f . As future work, we plan to explore using a bigger number of points DBW j for w > 3 to obtain a better approximation of complex 3Dworkpiece surfaces for the use in robotic automated fiber placement applications.