Global optimization of functional redundancy in a 6R robot for smoothing five-axis milling operations

Functional redundancy optimization is one of the most effective ways to improve robotic posture and performance. However, owing to the trade-off between smoothness of the robot's motion and other important aspects of performance, the optimization problem is too nonlinear to be solved in robotic milling processes. In this article, a derivation is shown that the robot's motion will be approximately smooth when the functional redundancies for each cutter location are smoothing variables. On this basis, a Legendre surrogate model is constructed to plan the overall optimal functional redundancy by remarkably reduced computation. Moreover, a trimmed sequential linear programming method is proposed to improve the solution further. From the outcomes of a 3D S-shaped grooves real milling experiment, the robot motion paths planned by the proposed methods achieve more stable tool-tip milling performance than that of comparable methods, while being several times or even 10 times faster in computational efficiency.


Introduction
The rapid development of modern manufacturing has made great demands on customization and intelligence in machining production (Ji and Wang 2019;Ji, Yin, and Wang 2019).To meet these requirements, robotics-based manufacturing has been increasingly applied to various machining tasks owing to its low cost, multi-functionality and large workplace applicability (Verl et al. 2019;Kim et al. 2019;Chen and Dong 2013;Pham and Ahn 2021).However, the robotic performance is still unstable, being limited by insufficient rigidity, poor accuracy and complex programming (Karim and Morris 2013).Functional redundancy optimization is one of the most effective ways to improve robotic machining performance and is widely applied in many five-axis machining operations such as robotic welding (Huo and Baron 2011;Zhou, Wang, and Gu 2022), drilling (Zhu et al. 2013;Wang et al. 2022) and milling (Zargarbashi, Khan, and Angeles "Posture Optimization in Robot" 2012;Xiao and Huan 2012;Zargarbashi, Khan, and Angeles "The Jacobian Condition Number" 2012;Léger and Angeles 2016;Yao et al. 2022).While functional redundancy refers to the rotation of the end effector along the machining direction for 6R industrial robots and has no direct impact on the machining target, its optimization is necessary and fruitful for extending closed-form solutions of robotic inverse kinematics (Shahabi and Ghariblu 2020;Peng et al. 2020;Lu et al. 2022).Since a robot's machining performance is closely related to its posture (Huo and Baron 2008;Bu et al. 2017;Lin, Zhao, and Ding 2017;Feng et al. 2021), the optimization of functional redundancy could improve the robotic machining performance.However, solving the optimization problem is susceptible to falling into local optima solutions in some continuous machining applications like robotic milling.This process involves the resolution of functional redundancy for a set of cutter locations and the trade-off between multiple machining performances, such as smoothness of the robot's motion, singularity and stiffness, which has need of a multi-objective optimization problem with many variables.
Many proposed methods are significant for solving the optimization problem in many robotic machining operations.W. Zhu et al. (2013) proposed a singularity performance index to evaluate the robotic motion dexterity at the joint displacement level to improve robotic drilling accuracy by optimizing the functional redundancy numerically.Yin et al. (2017) proposed a Cartesian compliance model to describe the stability of robotic drilling and chose a proper robotic drilling posture by repeatedly rotating the functional redundancy.Lin, Zhao, and Ding (2017) obtained a deformation map to determine the optimized robotic machining posture using proposed performance indexes like robot stiffness.Moreover, Jiao et al. (2018) proposed a combination of the singularity performance index and the robot static stiffness model to improve robot machining quality by planning the functional redundancy resolution.While these optimization methods involved the resolution of functional redundancy for only a single cutter location, the formulation of multiple performance criteria could still be a guide for other robotic machining tasks.
Apart from robotic drilling operations, other methods also play a critical role in optimizing functional redundancy for continuous machining tasks like robotic welding (Huo and Baron 2011).With this process, performance indexes such as robot motion smoothness are improved owing to the increased number of cutter locations.Shibata et al. (1997) introduced the weighted sum of squares of the joint acceleration as the global criterion, and the resolution of the functional redundancy using a genetic algorithm was significant in smoothing the robot's motion path.However, this method is not competitive in terms of computational cost and convergence.Dolgui andPashkevich (2009, 2004) suggested a graph-based method in the shortest path optimization scheme that has been successfully applied to smoothing a robot's motion path for high-speed robotic laser cutting.However, the graphbased method is still time-consuming with more cutter locations.Huo and Baron (2011) proposed that the self-adaptation weights could balance the performance indexes of joint limits and singularity avoidance.But the performance index smoothness of the robot's motion path was ignored, and the method has difficulty achieving robotic welding with high stability.
In robotic milling, the optimization problem is much more complex because, in the machining task, not only should the path be continuous, but also so should tool-tip contacts to the workpiece surface.This process increases the trade-off difficulty of having multiple performance indexes.Xiao and Huan (2012) advanced a focused algorithm for post-processing, solving functional redundancy for robotic motion planning, and achieving high performance.Other potential performance indexes can also be integrated in this method.However, the nonlinear search procedures are hard to characterize as highly accurate.Xiong, Ding, and L. Zhu (2019) introduced a robotic stiffness performance index, and the resolution of functional redundancy was used to improve machining stability and accuracy in the robotic milling processes.However, the smoothness performance index (Sun et al. 2016(Sun et al. , 2019;;Xu, Zhang, and Sun 2019) was ignored in robotic motion planning and the optimization problem was still nonlinear.Recently, Peng et al. (2020) developed a global optimization methodology to trade-off the smoothness of the robot's motion and other important performance indexes, such as singularity, stiffness and joint limit.Additionally, the Sequential Linear Programming (SLP) scheme has been successfully used to solve the functional redundancy accurately for high quality robotic milling operations.However, limited by the local iterative framework, the smoothness-oriented method critically relied on the initializations produced by the graph-based method (Dolgui and Pashkevich 2009), which was susceptible to locally optimal solutions.In summary, owing to the trade-off between smoothness of the robot's motion and other important performance indexes, solving the nonlinear optimization problem of functional redundancy for the robotic milling operation is challenging.While the smoothness-oriented method is powerful to solve this problem accurately using the SLP scheme, it is still susceptible to falling into locally optimal solutions.In this article, a smoothness problem is formulated so as globally to optimize the robotic motion path considered with multiple performance indexes.Then a derivation is shown that the robot's motion will be approximately smooth when the functional redundancies for each cutter location are smoothing variables, and a new smoothness performance index is proposed to simplify the optimization problem.Thus, a one-dimensional Legendre surrogate model (Fan, Huang, and Wang 2014) could be constructed to approximate the constrained optimization problem.Since fewer cutter locations and corresponding functional redundancies need only be uniformly sampled for constructing the Legendre surrogate model, the global optimization problem of functional redundancy for robotic milling processes can be solved efficiently.Moreover, owing to the reconstruction error of this model, a trimmed SLP method is proposed to locate these poor functionally redundant solutions dynamically and to program them linearly since the trimmed SLP reduces computational cost in each iteration.
Section 2 formulates the optimization of functional redundancy as the smoothness problem of a robot motion path constrained using multiple performance indexes for robotic milling operations.Section 3 analyses the correlation between the robot joint variables and the functionally redundant variables.On this basis, the Legendre surrogate model is constructed to approximate and solve the global optimization of functional redundancy fast.Section 4 proposes a trimmed SLP to improve the functional redundancies produced by the Legendre surrogate model.Section 5 uses several simulations and actual milling experiments to verify the effectiveness of the proposed methods.Section 6 summarizes the conclusions of this article and future work.

Problem formulation
The optimization problem of functional redundancy is formulated as the smoothness of a robot motion path constrained using multiple performance indexes based on the smoothness-oriented method (Peng et al. 2020).Given a series of discrete cutter locations CL = {cl i = (c i , z i ), i = 1, 2, . . ., N}, where c i ∈ R 3 and z i ∈ R 3 denote the position and the unit vector along the toolaxis direction in the robot base frame, respectively.Closed-form solutions of the robotic inverse kinematics are first expressed as the following: where η i ∈ [−π, π] denotes the functionally redundant variable of the robot end-effector; q i = {θ j i , j = 1, 2, . . ., 6} denotes the jth joint variable when the industry robot is milling at the ith cutter location; and θ j i denotes a vector of closed-form solutions.The model of robotic milling at the ith cutter location is shown as Figure 1. 1  For each specific η i and cl i , the number of θ j i 's is equal to 8 or 16 (Manseur and Doty 1989).Since the resolution of η i is discretely searched from the interval [−π, π] for each cutter location cl i , the number of robotic milling postures increases to infinity (Peng et al. 2020).Thus, the robotic milling posture could be flexibly planned for each cutter location, as shown in Figure 1.To offer high-quality robotic milling performance, the resolution of the functional redundancy is carefully determined by minimizing the following constrained formulation: where smooth = N−1 i=2 φ acc,i denotes the smoothness performance index of the robot motion path, and φ acc,i denotes the weighted squared sum of all joints' acceleration, expressed as follows: where ω j denotes the balanced weight and a 2 i,j denotes the jth joint's motion acceleration, which can be expressed as (Gautschi 2012): where t i = cl i+1 − cl i /f denotes the machining time with the commanded feed rate f.Moreover, Equation ( 2) is constrained by the performance indexes of joint limits and through singularity and stiffness to improve robotic milling performance.Among them, q min = {θ j min , j = 1, 2, . . ., 6} and q max = {θ j max , j = 1, 2, . . ., 6} denote the lower and upper bounds to the joint limits, respectively.κ(J(q i )) denotes the singularity performance index expressed by the condition number of the Jacobian matrix (Zhu et al. 2013): where is the normalization of the robot Jacobian matrix J(q i ) at the origin of the robot-end frame; L denotes the robotic characteristic length; • denotes the Frobenius norm.Since a smaller condition number κ(J(q i )) means the robot location is further from the singularities, a threshold δ 1 is preset to achieve robotic posture with good dexterity.Additionally, the stiffness performance index is computed to improve the robotic milling stability, expressed as follows (Guo, Dong, and Ke 2015): where det(C tt (q i )) denotes the determinant of the translational compliance matrix (C tt (q i )) at the tool end.Contrary to the singularity performance index, the larger is the condition number K stiff (q i ), the smaller is the deformation of the tool tip affected by the robotic milling force.For this reason, another preset threshold δ 2 is flexibly used to weight the maximum value K i s max and minimum value K i s min of the stiffness performance indexes for each cutter location.In this way, the optimization problem of functional redundancy has been formulated to improve the machining performance of the robotic milling operation.
While the optimization of functional redundancy for robotic milling is essentially a multi-object optimization problem with many variables, Equation (2) actually forms a constrained optimization formulation to achieve the resolution of functional redundancy for improving multiple performance indexes.However, the objective function involves six coupled joint variables at each cutter location, which is still too nonlinear to be optimized efficiently.The next section presents a simple and novel algorithm to address this problem.

Legendre-surrogate-model-based global optimization method
In order to simplify the optimization problem, a derivation is first shown demonstrting that the robot's motion will be approximately smooth when the functional redundancies for each cutter location are smoothing variables, and a new smoothness performance index is proposed for the robotic milling operations.On this basis, a more intuitive Legendre surrogate model is constructed to approximate the optimization problem, thus it is able to solve the overall optimal resolution of the functional redundancies fast and plan the robot's motion directly.

Analysis of the smoothness performance index for the robotic milling operation
In this section, a correlation between the functionally redundant variable and the robot joint variable is analysed to simplify the smoothness performance index for the robotic milling operation.Using the partial differential theory of the joint motion, the relationship between the sum of squares of acceleration φ acc,i and the functional redundancies η i for each cutter location cl i is expressed as (Peng et al. 2020): where ∂q i /∂η i denotes the partial differential between the joint variable and the functionally redundant variable for a specific cutter location, expressed as where J i Tcp (q i ) denotes the Jacobian matrix in robotics at the tool centrepoint, and is a 6 × 1 vector that means the functional redundancy affects only rotation of the tool tip along the machining direction.In order to derive the correlation of these variables for adjacent cutter locations, Equation ( 8) is calculated as a the definite integral: Given a specific Jacobian matrix J i Tcp (q i ), the joint variables q i are equivalent to: After that, each feedrate of the joint variables q i , i = 2, 3, . . ., N, and the machining time t i−1 are expressed as: Here, the joint variables of the robot motion path are assumed to be smoothed, and Jacobian matrices solved by the joint variables are similar for adjacent cutter locations.For this reason, Equation ( 11) is derived as approximately: As shown above, Equation ( 12) expresses the correlation between the functionally redundant variable and the robot joint variable in robotic milling processes when the robot joints' motion is assumed to be smoothed, and the resolution of the functional redundancy for each cutter location is also changed smoothly.Thus, a new smoothness performance index is obtained when robotic milling is achieved managing functional redundancies.In support of the derivation, the resolution of the functional redundancy for each cutting location is considered to be a smooth change .In fact, the joints' motion is approximately continuous and smooth in robotic milling processes.The reason is that the position and the machining directions of the adjacent cutter locations are similar to each other.Given the similar functional redundancies for adjacent cutter locations, the joint variables are similar according to robotic inverse kinematics, and result in the smooth performance of the robot motion path.

The construction of the Legendre surrogate model
As discussed in Section 3.1 regarding the smoothness performance index, the smoothness of the joints' motion is equivalent to the smooth functionally redundant variables for each cutter location.For this reason, the optimization problem of functional redundancy as shown in Equation ( 2) is simplified using the smooth Legendre surrogate model LG(x(c i )), revealing the following: where x(c i ) ∈ [−1, 1] denotes a one-dimensional variable mapped by the three-dimensional position c i of each cutter location.In the dimension reduction process, the Euclidean distance d i of adjacent cutting locations is expressed as Meanwhile, the Legendre surrogate model is expressed as follows (Fan, Huang, and Wang 2014): LG To solve the correction coefficient, the least squares method is combined with Equation ( 15), which is expressed as: where C sample = {c p ∈ CL, p = 1, . . ., M and M ≤ N} denotes the sampling positions within the cutter locations; η sample = {η p , p = 1, . . ., M} denotes the corresponding response values; β = {b 0 , . . ., b k d } denotes a vector of the correction coefficients; and LG(X(C sample )) denotes the basic matrix, which is expressed as follows: LG Then, the correction coefficient vector is solved by: From Equation ( 19) based on the specific sampling positions C sample and corresponding response values η sample , the functional redundancies for each cutter location are easily solved using the Legendre surrogate model.In addition to the smoothness performance index of the robot, the improvement of other performance indexes are considered to achieve a highly stable and accurate robot milling operation.For this problem, the solutions of the sampling positions and corresponding response values are developed in the next section.

The Legendre surrogate model to plan the overall optimal functional redundancy
To optimize the robot motion path globally, the sampling positions are uniformly sampled from the original cutter locations CL to produce C sample = {c p , p = 1, . . ., M and M ≤ N}.The floor(•) denotes the rounded element to the nearest integer toward the direction of negative infinity.Based on this, the functional redundancies are solved as response values for each sampling position.Since the multiple performance indexes in the constraints of Equation ( 13) are too loose to solve a specific functional redundancy, this equation is rewritten as follows: The stiffness of the constraint performance index in Equation ( 20) is the maximal value, while the singularity performance index is constrained within a good interval Owing to the relatively robotic low stiffness compared to the computer numerical control machine, used to automate, control and monitor the movements of a machine.Moreover, deviation of the positioning accuracy, poor machining quality and inferior production efficiency could be induced by cutting forces (Guo, Dong, and Ke 2015;Zhang et al. 2005;Olsson et al. 2010;Liao et al. 2021).At present, search space discretization (Dolgui and Pashkevich 2009) is introduced to solve the response values.The parameter π] is also uniformly sampled with the step b = 2π/m k , where m k denotes a preset positive integer.Then, the m k + 1 sampling nodes within the feasible interval [−π, π] are expressed as η t p = −π + t • b, t = 0, 1, . . ., m k .For each sampling position c p , the sampled nodes η t p are tested by the constraints of Equation ( 20).Moreover, the best sampled node η t best is regarded as the response value η p for each sampled cutter location when the feasible robotic posture performance has good stiffness and is far away from singularities.To this end, by combining the sampling positions C sample and the corresponding response values η sample = {η p , p = 1, . . ., M}, the Legendre surrogate model can be constructed to plan the overall optimal functional redundancy and the robot motion path.15), ( 17) and ( 19) with x(c p ) and η sample .
4 Mapping the position c i of each cutter location cl i as x(c i ) by Equation ( 14). 5 Computing the globally optimized functional redundancy η opt by Equation ( 13) with the constructed Legendre surrogate model LG(x(c i )) and x(c i ).6 Computing the robot joint variable q opt by the inverse kinematics and Equation ( 21) with η opt and cl i As shown in of Algorithm 1, the Legendre-surrogate-model-based global optimization method presents the details to solve the functional redundancy for the robotic milling process fast.The principle is simple and effective, as shown in Figure 3 at https://doi.org/10.6084/m9.figshare.21444312,which includes: (1) developing a set of cutter locations; (2) constructing a one-dimensional Legendre surrogate model; and (3) solving optimized functional redundancy η opt = {η opt i , i = 1, 2, . . ., N} for each cutter location.Lines 1-4 of Algorithm 1 show that search space discretization needs to be used only for sampled cutter locations, aimed at solving optimized functional redundancy (response values) for constructing Legendre surrogate model.Thus, the computational cost of this process is lower even if the number of cutter locations is increased for the robotic milling operation.After that, given a series of discrete cutter locations, Lines 5-7 of Algorithm 1 show that the functional redundancies and joint variables can be solved directly to plan the robot's motion path efficiently.It is noted that the closed-form solutions of the inverse kinematics solved by using Equation (1) are still not unique (Manseur and Doty 1989), which should be arranged as follows (Peng et al. 2020;Dolgui and Pashkevich 2009): where q opt = {θ opt i,j , i = 1, 2, . . ., N, j = 1, 2, . . ., 6} denotes the optimized joint angles of the robot.For this problem, dynamic programming (Dolgui and Pashkevich 2009) is one of the most effective methods for solving Equation (21).

Trimmed sequential linearization programming method
In this section, while the resolution of the functional redundancies has been overall optimized by the Legendre surrogate model in the robotic milling processes, its results can still be improved to achieve better solutions owing to the reconstruction errors of this model.Thus, this article linearizes the optimization problem of functional redundancy and proposes a trimmed SLP scheme to improve the solutions produced by the Legendre surrogate model efficiently.

Linear formulation corresponding to the functional redundancy optimization
Given a series of functional redundancies produced by the Legendre surrogate model, it can be further improved using a linear programming scheme.Based on the smoothness-oriented method (Peng et al. 2020), the constrained optimization problem in Equation ( 2) is firstly linearized as follows: where ml denotes the imposed move limit of the functionally redundant variable, is the feasible subinterval to the functional redundancy variable η i ∈ [−π, π] for each cutter location.Based on the search space discretization method discussed in Section 3.3, the η i is uniformly sampled as η t i = −π + t • b, t = 0, 1, . . ., m k with the step b = 2π/m k .For each cutter location, the sampled nodes η t i are tested by the constraints of Equation (2).After that, the feasible sampling nodes η k f i ⊆ η t i are regarded as the lower and upper bounds of the feasible subintervals . This is obviously a standard linear programming problem that can be solved by the conventional Simplex algorithm.However, the computation costs will increase with an increase in the number of functionally redundant variables.For this problem, it can be seen that most of the functional redundancies produced by the Legendre surrogate model were good.Only a few of them must be improved, which helps to reduce the computational cost.For this purpose, a trimmed sequential linear programming scheme is presented in the next section to solve this problem efficiently.

The trimmed SLP for the robot motion planning
Following the trimmed operation developed by Chetverikov, Stepanov, and Krsek (2005), it is possible to trim the poor functional redundancies produced by the Legendre surrogate model.The weighted sum of squares of the joints' acceleration φ acc,i , i = 2, . . ., N − 1, is sorted in increasing order.Then, the least trimmed value N po is given to select the rear N − N po − 2 elements from the sorted φ acc,i,i * , i * = 1, 2, . . ., N(φ acc,i,1 ≤ φ acc,i,2 ≤ . . .≤ φ acc,i,N−3 ≤ φ acc,i,N−2 ).While selecting as small a φ acc,i as possible, the optimal parameter N po is automatically set by minimizing the energy function as follows: where λ ≥ 0 is a preset parameter, and e(N po , φ acc,i,i * ) denotes as the trimmed mean squared error function: where N po i=1 φ acc,i,i * denotes the summing of the first N po elements for the sorted φ acc,i,i * .As discussed in Section 4.2 regarding Chetverikov, Stepanov, and Krsek (2005), ψ(N po , φ acc,i,i * ) is a smooth and convex function that has a distinct single minimum for which the modified golden section search algorithm is used to solve the optimal parameter N po .In this study, the goal is to minimize Equation ( 21) for the resolution of the functional redundancy.The key observation is that some resolutions of the functional redundancies produced by the Legendre surrogate model are better.Thus trimming these functionally redundant variables reduces the computational costs and improves the sequential linear programming.For this reason, the modification is the initial φ acc,i , i = 2, . . ., N − 1,inthefirstiteration.Theabovetrimmedoperationsortsandselectsthepoor φ acc,i,i * (i * = N po , . . ., N − 2) to formulate the sequential linear programming: where i(i * ) denote an array index of cutter locations that correspond to a higher value of acceleration φ v acc,l oo .Functionally redundant variables η i(i * )−1 , η i(i * ) , η i(i * )+1 are repositioned to wait for further improvement.On the contrary, most of the functionally redundant variables corresponding to lower values of acceleration φ acc,i,i * (i * < N po ) are set to zero and eliminated through the trimming operation.On this basis, Equation ( 25) is a standard sequential linear programming problem, which can be solved using the traditional Simplex method.As shown in Algorithm 2, the functional redundancy produced by the Legendre surrogate model is given as an initial solution.With iterations, Algorithm 2 alternates between updating the functional redundancies under the current optimization problem and trimming the solutions to formulate the optimization problem of trimmed sequential linear programming until convergence is achieved, as shown in Figure 4 at https://doi.org/10.6084/m9.figshare.21444312.
Similar to the smoothness-oriented method (Peng et al. 2020), a strategy is applied to shrink the imposed move limit ml dynamically to achieve accuracy of the sequential linear programming.While Equation ( 7) performances the gradient to improve the computational efficiency of the SLP for this procedure, the move limit makes it valid.The smaller is the move limit, the more accurate is the gradient to be used for solving Equation ( 25), but the convergence speed slows down.The dynamical Algorithm 2: Functional redundancy optimization based on the trimmed Simplex method in the robotic milling process Input: Cutter locations CL = {cll i = (c i , z i )}, i = 1, 2, . . ., N, the initial functional redundancies η opt = {η opt i , i = 1, 2, . . ., N} produced by the Legendre surrogate model of Algorithm 1, the commanded feed rate f , the preset parameter λ, the number of the sampling nodes m k + 1, threshold values ε 1 , ε 2 and a maximum number of iterations IterMax.Output: Optimized solution η * = {η * i , i = 1, 2, . . ., N} and corresponding joint variables q * = {θ * i,j , i = 1, 2, . . ., N, j = 1, 2, . . ., 6}.

else 8
Quit the loop.9 end 10 end 11 end move limit offers a trade-off between solution quality and computation time that is also used as one of the stop criteria.If the current objective function value is larger than the previous one, the move limit is decreased, otherwise it is increased, but to no larger than the initial preset value.For convergence, when the move limit is smaller than the preset threshold ε 1 , another threshold ε 2 is larger than the magnitude between the current objective function value and previous ones, or the iteration numbers are larger than the maximum allowable one, the optimization of the functional redundancy is finished.Algorithm 2 presents the details.

The experimental results
In this section, the performance of the following proposed methods-(Version 1, V1), (Version 1+Version 2, V1+V2), graph-based (Dolgui and Pashkevich 2009) and smoothnessoriented (Peng et al. 2020)-are compared through simulation cases and real milling experiments.The proposed V1 is the Legendre-surrogate-model-based global optimization method given in Section 3, and the V1+V2 method denotes the trimmed SLP shown in Section 4 to improve the solutions of the proposed method V1.In addition, the KUKA kr210-R2700 industrial robot is equipped in the robotic machining system.The joint limits of the robot are q min = [−37π/36, −5π/18, −7π/6, −35π/18, −25π/36, −35π/18], q max = [37π/36, 17π/36, 13π/36, 35π/18, 25π/36, 35π/18].Following the stiffness identification experiment developed by Dumas et al. (2011) and Cordes and Hintze (2017), the robotic characteristic length is solved as L = 708 mm, and the joint stiffness vector of the robot is defined as k θ = [1.56, 1.63, 3.45, 1.23, 0.94, 0.18] × 10 9 N mm/rad.Additionally, the highest order of the Legendre polynomials is set to six, and the number of sampling points is set to 60 to construct the Legendre surrogate model.For the trimmed SLP and the smoothness-oriented method, the constrained weight of the smoothness, the singularity and the stiffness index are set to ω i = 1 (i = 1, 2, . . ., 6), δ 1 = 36 and δ 2 = 0.9, respectively.The step size of the discrete search is set to 2π/M, M = 180; the thresholds for convergence are set to ε 1 = 0.0001π/180, ε 2 = 0.001.The preset parameter for the trimming operation is set to λ = 2; the maximum number of iterations is IterMax = 150, and the imposed move limit ml = 0.5 × 2π/M.After these preparations, the comparison methods were implemented on a standard computer with the following specifications: a 3.4 GHz Intel Core tm i7-6700 CPU with 12 GB RAM through MATLAB code.

Results of S-shaped trajectory planning
The commercial CAM software, UG NX, is used for simulation and to generate 3D S-shaped grooves in the world frame, as shown in Figure 5 at https://doi.org/10.6084/m9.figshare.21444312.For the tool-end trajectories of robot milling, the 3D S-shaped grooves produced by a standard cylinder have different intercepts of D = 60 mm, with the corresponding number of cutter locations N = 300.Moreover, the tool-axis directions for each cutter location are defined by the unit normal vectors of the corresponding cutting positions on the cylinder surface; the transform from the robot's base frame and the world frame is set to T W B = (2000, −200, 600, cos(π/8), 0, 0, sin(π/8)).Then, the transformation relationship between the flange frame and the milling tool frame Tcp is set to T Tcp 6 = (300, 0, 100, cos(π/4), 0, 0, sin(π/4)), where the first three elements are the position (in millimetres), and the other four elements denote the unit quaternion.Furthermore, the robotic kinematic model is expressed as follows (Peng et al. 2020): (26) where T x (π) denotes a homogeneous matrix with a rotation angle of π around the x-axis; T Tcp 6 denotes a homogeneous matrix of the robotic forward kinematics, calculated from the DHparameters shown in Table 1 at https://doi.org/10.6084/m9.figshare.21444312;f −1 (•) denotes the achieving of closed-form solutions of the robotic inverse kinematics; T z (η i ) denotes a homogeneous matrix with a rotation angle η i around the z-axis.The corresponding local coordinate system is defined by the cutter locations CL = {cl i = (c i , z i )}, i = 1, 2, . . ., N, where the z-axis is defined as z 0 i = z i ; the y-axis is defined as Combining with these kinematic parameters, the robot feed rate command is set to f = 600 mm/min.Then, comparative methods are used to plan the motion path for the robotic milling processes.As shown in Figure 6 at https://doi.org/10.6084/m9.figshare.21444312,it can be seen that the joint angles with the graph-based method have fluctuations locally, and the smoothnessoriented method smooths these fluctuations to obtain better results.However, by comparing the proposed method V1 and the proposed method V1+V2, the joint angles planned using the smoothnessoriented method were not flattering enough, resulting in worse smoothness of the robot motion.As shown in Figure 7 at https://doi.org/10.6084/m9.figshare.21444312,the accelerations of joint 4 and the weighted sum of squares of all joints' acceleration are used to evaluate the smoothness performance index of the robot motion.The maximum (average) accelerations from the proposed methods are 6.12e-02 (1.52e-02) rad/s 2 (V1) and 6.46e-02 (1.46e-02) rad/s 2 (V1+V2).Compared to the graph-based method 3.24 (0.34) rad/s 2 and the smoothness-oriented method 0.12 (5.26e-02) rad/s 2 , the results of the proposed methods decreased by approximately 98% (95%) and 49% (71%), respectively.Meanwhile, the maximum (average) smoothness indexes are 7.00e-03 (7.23e-04) rad/s 2 and 6.30e-03 (6.42e-04) rad/s 2 for the proposed methods V1 and V1+V2, which decrease by approximately 99% (99%) rad/s 2 and 61% (90%) rad/s 2 compared with the graph-based method (16.29 (1.19) rad/s 2 ) or the smoothness-oriented method (1.81e-02 (7.40e-03) rad/s 2 ).
The results showed that the proposed methods V1 and V1+V2 performed better than the other methods for smoothing the robot motion.The reason is that the discrete search process of the graphbased method constrained by all-round performance indexes is susceptible to trap in local minimum, resulting in poor solutions for the functional redundancy.Since the smoothness-oriented method relies on the solution of the graph-based method, its solutions also easily plunge into the local optimum together.In contrast, the proposed method V1 improves all-round robot performance indexes, and allows the simpler Legendre surrogate model to smooth the functional redundancy effectively for the robotic milling operation, as shown in Figure 8 at https://doi.org/10.6084/m9.figshare.21444312.It can be seen that the overall stiffness performance of the corresponding joint configuration by the proposed methods are also better than the other methods compared.While the singularity condition numbers corresponding to the proposed methods are slightly larger, the maximum number is 2.75, denoting that the optimized robot motion path is far from being a singularity.
In terms of efficiency, the running times of the proposed methods V1 and V1+V2 are 7.02 s and 10.37 s, respectively, which are 3-5 times faster than the graph-based method (33.46 s) or the smoothness-oriented method (37.38 s).The computational cost of the proposed method V1 is partly the reason for the optimization process of functional redundancy for sampled cutter locations.However, the graph-based method discretely searches for reasonable functional redundancy for each cutter location.Since the proposed method V1+V2 and the smoothness-oriented method aim to improve the solutions produced by the proposed method V1 and the graph-based method, respectively, their computational cost depends on these global methods.

Real applications
As discussed previously in Section 5.1, the comparative methods improve the robotic machining stiffness, and keep the robot posture away from the singularity.The performances of the proposed methods are better than those of the graph-based method and the smoothness-oriented method regarding the smoothing of the robot motion path.In addition, a machine-vision-guided robotic milling processing platform is built to verify the effectiveness of the proposed methods, as shown in Figure 9 at https://doi.org/10.6084/m9.figshare.21444312.Herein, the cylindrical workpiece material is ductile iron QT500-7; the weight is appromaxitely 90 kg; the electric spindle model is DJ110-9.0-12.0KW-ISO30.The command feed rate increases to f = 900 mm/min; the speed of the spindle is set to 12,000 rev/min; and the end of the tool is equipped with a tapered tungsten steel cutter to mill the workpiece.For the location problem, the transformation from the flange to the tool end is first calibrated as T Tcp 6 = (357.86,0.09, 96.02, cos(π/4), 0, 0, sin(π/4)).Then, a PS801-E1 depth camera is used to sensor the workpiece to achieve a 3D point cloud in the robot base frame.Moreover, the S-shaped groove (D = 120 mm, N = 500) is located in the same coordinate system.By combining human-machine interaction with an Iterative Closest Point algorithm (Besl and Mckay 1992), the cutter locations on the S-shaped groove are transformed from the world coordinate system to the robot base frame using the left transformation T best = (1674.43, 22.32, 638.33, −0.688, 0.0118, 0.0047, −0.725).
As shown in Figure 10 at https://doi.org/10.6084/m9.figshare.21444312,the joint angles are planned for the robotic milling operation on the 3D S-shaped groove based on the outcomes of different comparative methods.The joint angles with the graph-based method are shaped zigzag fluctuations locally, while the smoothness-oriented method also failed to improve these poor solutions significantly.On the contrary, the results of the proposed methods are more smooth than those of the other methods.The maximum (average) smoothness indexes of the joint motion are 0.25 (1.14e-02) rad/s 2 and 9.70e-02 (0.87e-02) rad/s 2 for the proposed method V1 and the proposed method V1+V2, respectively.Compared with the graph-based method, 629.07 (103.41)rad/s 2 , and the smoothness-oriented method, 317.47 (68.59) rad/s 2 , the reductions in the smoothness index are more than 99%.The minimum (average) stiffness indexes and maximum (average) singularity conditions of the proposed methods are about 703.71 (790.81)N/mm and 2.59 (2.40), respectively, which are better than the results-693.12(779.74)N/mm and 2.82 (2.47)-of other methods compared, as shown in Figure 11 at https://doi.org/10.6084/m9.figshare.21444312.Moreover, the running times of the proposed methods V1 and V1+V2 are 6.97 and 14.65 s, respectively, less than the graph-based method (53.68 s) and the smoothness-oriented method (58.50 s).From the simulation results, the proposed methods are more efficient than the other methods compared for optimizing the robot motion path.
Following the simulation results, the KUKA kr210-R2700 robot was used to mill the cylindrical cast iron: the diameter of the bottom surface was 150 mm; the height was 500 mm; and the milling depth is 0.5 mm.The machined S-shaped grooves using the comparative methods are shown in Figure 12 at https://doi.org/10.6084/m9.figshare.21444312.The results show that the surface quality of machined S-shaped grooves with the proposed methods are better than those of the graph-based method and the smoothness-oriented method.The surface quality of the machined S-shaped grooves with the graph-based method and the smoothness-oriented method have led to unfavourable wavy marks.The proposed methods are much better, considering the actual milling efficiency.The machining time of the proposed methods are approximately 37 s.The machining times of the graph-based method and the smoothness-oriented method are approximately 45 s.The machining paths planned by the proposed methods, compared with the other methods, provides efficient smoothness of the robot joint motion and the milling task.
While the performances of the proposed methods V1 and V1+V2 are pretty similar or little different for this real milling case, the proposed method V1+V2 is still contributive and valuable.In this article, the proposed method V1 (the Legendre-surrogate-model-based method) aimed to solve the overall optimal resolution of functional redundancy globally for robotic milling processes.Additionally, the proposed method V1+V2 is presented to add the proposed method V2 (trimmed sequential linear programming) as a local iterative solver to improve the solution further produced by the proposed method V1.In this process, the better is the solution of the proposed method V1, the smaller will be the improvement of the proposed method V2.Thus, the performances of the proposed methods V1 and V1+V2 should be pretty similar when the solution of the proposed method V1 is pretty good.However, the proposed method V1+V2, which aimed to improve the solution further produced by the proposed method V1, is more robust for other complex machining tasks.Moreover, the proposed V2 linearizes the highly nonlinear optimization problem and solves it by using trimmed sequential linear programming.This process makes the gradient optimization effective and makes it possible to achieve a solution of high accuracy.

Conclusion
Optimization methods are proposed to smooth a robot motion path constrained by multiple performance indexes in robotic milling processes.A global optimization method is first proposed to plan the robot motion path based on the Legendre surrogate model, and then a derivation formula is provided to show the correlation between the joint variables and the functionally redundant variables.To improve the solutions produced by the Legendre surrogate model, sequential linear programming is developed to formulate a constrained optimization problem of the functional redundancy that can be solved by using the trimmed SLP method based on an extra energy function.The optimization methods of functional redundancy achieve a better robot motion path than those of the graph-based method and the smoothness-oriented method.In a real milling experiment, the proposed methods provide fast and smooth solutions for the robot motion path.
The proposed methods can be applied to other continuous robot machining tasks integrated with different constraints, such as high-speed robotic welding and spraying.Moreover, the proposed method V2 is a new and efficient linear programming scheme that, significantly, can also be applied in other fields, such as portfolio investment, transportation and production organization.After that, a few directions can be used to improve the proposed method.Future studies should plan the motion path of the worktable or the position of the robot, which are powerful ways to extend the feasible robot motion path and machining performance.

Note
1. Figures 1-12 and Table 1 relating to this article are available free of charge on https://figshare.com at https://doi.org/10.6084/m9.figshare.21444312.Brief summary details: Figure 1 shows different robotic machining postures of the same cutter location.Figure 2 shows Legendre polynomials of the first six orders.Figures 3  and 4 show schematic diagrams of Algorithms 1 and 2, respectively.Figures 5-8 show the simulation results of robotic paths planned by four comparison methods on a virtual 3D S-shaped groove, and Table 1 is mentioned in the caption to Figure 5. Figure 9 shows the experimental platform of robot milling guided by machine vision.
Figures 10 and 11 show the simulation results of robotic paths planned by different comparison methods on a real 3D S-shaped groove.According to these results, Figure 12 shows the real milling results of the robot on the 3D S-shaped groove.
) where b k k (k d = 0, 1, . . ., N d ) denotes the correction coefficient of the Legendre surrogate model with the highest order of N d , and LG k d (x(c i )) denotes the Legendre polynomial of order k d as shown in Figure 2 at https://doi.org/10.6084/m9.figshare.21444312,and expressed as follows: