Optimal Robot Motion for Accelerated Life Testing of a 6-DoF Industrial Robot

In order to verify the reliability of drive components for industrial robots, component-level life tests must be accompanied by a system-level life test using actual robots in which predefined robot motions are repeated throughout the test. To properly verify the durability of drive components through a system-level life test, it is important to design test modes so that the required test time is the same for all joint drive components of the robot, and it is necessary to design test modes with a high acceleration factor so as to shorten the required test time as much as possible. To solve this problem, the present research proposes a method for designing robot motions that makes the accelerated life test time for all the drive components of the robot equal. In particular, we solve a dynamic based motion optimization problem for an industrial 6-DoF (degrees-of-freedom) robot that minimizes the AM-GM (arithmetic mean to geometric mean) ratio of the acceleration factors of each joint. The results show that C2-continuous test modes with the same acceleration factor, which is inversely proportional to the cycle time of the robot motion, can be derived.


Introduction
Industrial 6-DoF (degrees-of-freedom) robots are widely used for the automation of various manufacturing processes, such as welding, assembly, painting, inspection, etc. As the demand for industrial 6-DoF robots continues to increase, robot manufacturers are correspondingly researching the development and performance improvement of robots to meet the requirements of various applications.
It is imperative to ensure the performance and reliability of robot drive components including servo motors and gearboxes in improving the accuracy and responsiveness of industrial robots. In particular, a high-precision gearbox with little backlash and high reduction ratio is a key component which determines the accuracy and performance of a robot, and since the weight of the gearbox is directly related to the payload of robots, a small size, light weight, and high torque density design is required. Therefore, high-precision gearboxes, such as harmonic drives [1] and planocentric type 2-stage precision gearboxes [2], are applied to many industrial robots requiring high accuracy and responsiveness. High-precision drive components occupy the highest proportion in the manufacturing cost of robots, and it is essential to ensure stable performance and reliability because they can cause downtime and accidents of the automated process in the event of a failure.
In the field of reliability application, the most commonly used test methods for reliability qualification are zero-failure acceptance tests since they require fewer test samples and less test time compared to other test methods that guarantee the same reliability with a given confidence level [3]. To verify the reliability of drive components for industrial multi-joint robots, component-level life tests must be accompanied by a system-level life test in which the drive components are mounted on an actual robot, and predefined test modes are repeated to verify the durability.
A zero-failure acceptance test concept can be applied to calculate the test time for the guaranteed life of each joint part in system-level life test, and an acceleration test is designed by comparing the test load with the load in guaranteed conditions. The load applied to the drive components of each joint of the robot will depend on the payload and the motion profile. This can also result in differences in the degree of damage to the drive components, which results in a difference in the required test time.
As a consequence, the overall zero failure test time is determined by the longest test time among the required test times of each joint of the robot. If the variation of the required test time for each joint is large, the parts tested in excess of the required time have higher probability of failure, necessitating the repair or replacement of damaged components, and thus requiring more time and cost for the test. There is a method to stop the motion of the joint after the test time has been completed, but it can be cumbersome for the tester such as monitoring the time when the test time for each axis is completed and operating in a new test mode, etc. especially when the test scale is large.
Therefore, one of the ideas to verify the durable life of a robot drive component through system-level accelerated life testing is to design a test mode with the same zero failure test time for all the joint drive components. At the same time, it is necessary to design the test with a high acceleration factor to reduce the test time required to verify the durable life. For this purpose, the industry has been using robot simulation software to design test modes based on experience and intuition of engineers and trial and error. Therefore, it took a lot of time to derive robot motions that match the test times for the components of the individual joints equally, and there was a limit in the accuracy.
Motivated by the above observation, the goal of this study is to develop a numerical method for designing robot motions for system-level life test of drive components of robot.
In order to achieve the research goal, this study drew the following results: (1) A dynamic based optimization problem was formulated in which motion profiles defined over a fixed cycle time are used as optimization variables and the AM-GM ratio of the acceleration factors of each joint is used as an objective function. (2) For a given initial and final pose of a robot, a global minimum solution that that makes the accelerated life test times for the drive components of all joints of a robot equal was derived. (3) Using the kinematics and dynamics parameters of an actual 6-DoF industrial robot, a global optimal solution was obtained for various cycle times and it was confirmed that the acceleration factor is inversely proportional to the cycle time.
From the above results, it was possible to design robot motions which have the following useful properties for system-level life test: equal test time for the drive components of all joints of a robot, smooth joint displacement so that the required joint torque is continuous, and high acceleration factor to shorten the required test time.
Robotic motion optimizations for various robots have been studied to date [4][5][6][7][8][9][10], but most have focused on optimizing efficiency and performance, such as minimizing energy consumption or shortening cycle time to perform a particular task. What distinguishes the present study from past works relating the robot motion optimization is that the present study is a result of studies on motion optimization to ensure that the load applied to each drive component is evenly distributed.
In refs. [4][5][6][7][8], motion optimization problems are established with robot dynamics based on Lie group. A recursive algorithm to obtain the analytic gradient that can be used to perform motion optimization is presented, and the displacement of the robot joint was parameterized with a B-spline curve to convert the problem into a finite-dimensional constrained optimization problem. In Martin et al. [6] and Bobrow et al. [7], minimum effort optimal motions for various open chain manipulators are examined, In Wang el al. [8], the result of applying motion optimization for maximizing payload to the industrial robot PUMA is presented. The result of applying motion optimization technique to the hydraulic excavator and the bipedal humanoid robot are given in [9] and [10], respectively. The efficient generation of the natural robot motion based on dynamics-based optimization and principal component analysis is presented in [11].
Studies on the reliability and lifetime of robots have been conducted mainly on a component-by-component basis, and there have been few studies to systematically design the life test at system-level. In this paper, it was attempted to design robot motions for the accelerated life test of drive components at system level by applying a robot motion optimization technique.
Studies on the reliability and lifetime of precision gearboxes used in robots are presented in [12][13][14] and work on servo motors are summarized in [15]. In ref. [16], an accelerated life test was performed on individual components to calculate the reliability of the entire system from the reliability of each component of the robot, and the average lifetime of the system was estimated from the results. In ref. [17], a study on motion optimization to increase the service life of the robot gearbox was conducted.
This paper is configured as follows. In Section 2, modeling for the dynamic analysis of industrial 6-DoF is presented. Section 3 introduces a concept of a zero failure life test and an acceleration factor in relation to an accelerated life test for drive components of an industrial 6-DoF robot. Section 4 establishes a motion optimization problem that can verify the lifetime of all the drive components within the zero failure test time by matching the acceleration factor of each robot joint equally. Section 5 shows the results of performing the motion profile optimization to the actual 6-DoF industrial robot. Section 6 summarizes the main results of the study as a conclusion. Figure 1 represents a 6-DoF industrial robot that is the subject of this study. The robot link is numbered from 1 to 6 and the body-fixed coordinate frames are attached to each link as in the figure. The coordinate frame 0 is a space fixed frame of reference and coordinate frame 7 is attached to an object fixed to the link 6.

Dynamic Modeling of the Robot
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 16 optimization is presented, and the displacement of the robot joint was parameterized with a B-spline curve to convert the problem into a finite-dimensional constrained optimization problem. In Martin et al. [6] and Bobrow et al. [7], minimum effort optimal motions for various open chain manipulators are examined, In Wang el al. [8], the result of applying motion optimization for maximizing payload to the industrial robot PUMA is presented. The result of applying motion optimization technique to the hydraulic excavator and the bipedal humanoid robot are given in [9] and [10], respectively. The efficient generation of the natural robot motion based on dynamics-based optimization and principal component analysis is presented in [11]. Studies on the reliability and lifetime of robots have been conducted mainly on a component-by-component basis, and there have been few studies to systematically design the life test at system-level. In this paper, it was attempted to design robot motions for the accelerated life test of drive components at system level by applying a robot motion optimization technique.
Studies on the reliability and lifetime of precision gearboxes used in robots are presented in [12][13][14] and work on servo motors are summarized in [15]. In ref. [16], an accelerated life test was performed on individual components to calculate the reliability of the entire system from the reliability of each component of the robot, and the average lifetime of the system was estimated from the results. In ref. [17], a study on motion optimization to increase the service life of the robot gearbox was conducted.
This paper is configured as follows. In Section 2, modeling for the dynamic analysis of industrial 6-DoF is presented. Section 3 introduces a concept of a zero failure life test and an acceleration factor in relation to an accelerated life test for drive components of an industrial 6-DoF robot. Section 4 establishes a motion optimization problem that can verify the lifetime of all the drive components within the zero failure test time by matching the acceleration factor of each robot joint equally. Section 5 shows the results of performing the motion profile optimization to the actual 6-DoF industrial robot. Section 6 summarizes the main results of the study as a conclusion. Figure 1 represents a 6-DoF industrial robot that is the subject of this study. The robot link is numbered from 1 to 6 and the body-fixed coordinate frames are attached to each link as in the figure. The coordinate frame 0 is a space fixed frame of reference and coordinate frame 7 is attached to an object fixed to the link 6.   generalized velocity of the space fixed frame and the generalized forces exerted on the end-effector as in Equation (1). The load torque by gravity is reflected in the joint torque by initializing the generalized acceleration of the space fixed frame to the acceleration due to gravity.

Dynamic Modeling of the Robot
In the outward recursion phase, the coordinate transformation matrices between the adjacent links T i−1,i are calculated for the joint angle as Equation (2), where M i are the coordinate transformation matrices in home position. Then the generalized velocity V i and their time derivatives . V i are calculated recursively in the order from the base frame to the end effector as Equations (3) and (4).
In the inward recursion phase, the generalized forces F i for each link are obtained recursively by solving the Equation (5) in the order from the end effector to the base frame. The joint torque is then calculated from the joint motion direction component of the generalized force as Equation (6).
The derivation of the algorithm including the mathematical definition of the adjoint operators Ad T and ad V can be found in the original paper [18,19]. Initialization: Outward recursion (from i = 1 to i = 6): .

Accelerated Life Testing of Robot
Among the drive components of the high-precision robot dealt with in this study, the high-precision gearbox is the most prone to deterioration due to friction and wear, and thus can be seen as a component that has the greatest influence on the durability of the robot. The life of the high-precision gearbox is determined by the life of the bearing supporting the input-side rotary shaft. Therefore, the zero failure test time T n to verify the durability of the gearbox is derived from Equation (7) [3].
where B 100p is the guaranteed lifetime according to the conditions of use, CL is the required confidence level, n is the number of samples, p is the cumulative probability of failure, and β is the shape parameter. The guaranteed lifetime for the rated use conditions of the gearbox is usually given by component-level life test results from the manufacturer. The precision gearbox used in a variable load environment can then be guaranteed in service life only when used within the permissible torque presented by the manufacturer, which is related to the service life of the bearings of the precision gearbox components.
Ball bearings (e.g., in harmonic gearboxes) or roller bearings (e.g., in planocentric type 2-stage precision gearboxes) are used in precision gearboxes depending on the load class and characteristics [1,2]. Bearings have a well-established theory and wealth of verifications, and the lifetimes of precision reducers are also calculated according to this theory.

of 16
First, from the load torque τ and speed ω of the gearbox over a fixed cycle time of t f , equivalent load torque τ eq and equivalent speed ω eq are defined as Equations (8) and (9). Here, the index α is determined by the type of bearings used: α = 3 when ball bearings are used and α = 10/3 when roller bearings are used as input rotary parts [20].
Then, the acceleration factor AF of the motion profile having the equivalent load torque τ eq and the equivalent speed ω eq based on the rated torque τ rate and the rated speed ω rate of the robot's precision gearbox can be defined as the Equation (10).
Consequently, according to the definition of zero failure test time T n and the acceleration factor, the accelerated life test time T accel can be calculated as in Equation (11) [21].
The purpose of the present study is to design a test mode in which the system-level life test of precision gearboxes for robots is completed in the shortest time possible, and the durability of all gearboxes can be verified in one test without replacing samples during the life test.
In other words, the goal is to design a motion profile that allows the accelerated life test time T accel required by the precision gearbox mounted to each joint of the robot to be the same, and makes T accel as small as possible.

Formulation of Optimization Problem
Without loss of generality, and under the assumption that the zero failure life test time T n required for the gearboxes used in each joint of the robot are the same, all acceleration factors must be the same to ensure that the acceleration life test time T accel is the same for each joint.
All of the acceleration factors of the drive components derived from Equations (8)-(10) for each joint of the robot have non-negative values, so an inequality as shown in Equation (12) holds for their arithmetic and geometric means. Here, the degrees of freedom of robot is designated by N to make the expression more general. The necessary and sufficient condition to make them equal is that acceleration factors for all drive components have the same value.
Therefore, in this study, a function J, which represents the ratio of the arithmetic mean of the acceleration factors to their geometric mean, is defined as Equation (13), and its logarithmic value as in Equation (14) is defined as the cost function of the optimization problem. This ensures that the acceleration factor of each drive component converges to the same value when the global minimum is reached, where the cost function is zero as in Equation (15). log As a result, the problem of obtaining a motion profile in which all acceleration factor are the same during the predefined cycle time t f can be defined as a constrained optimization problem that minimizes the logarithm of AM-GM ratio (Equation (16)) while satisfying constraints on torque τ i (t), speed . q i (t), and displacement q i (t) of the robot drive components as in Equations (18)-(22). Here, the Equation (17) is the dynamic equation which is solved from the Equations (1)-(6) and τ max i , . q max i is limit on joint torque and speed, respectively. The lower and upper bound of each joint displacement is designated by q min i and q max i , and the initial and final value of the joint angles during the specific cycle time t f is q 0 and q f . .
Since the motion profile that causes the acceleration factor of the gearboxes of each joint to be the same is not unique the above optimization problem may contain multiple global minimum. Therefore, the solution of the above constrained optimization problem can converge to different global optimal solutions depending on the initial setting values of the optimization variable. Moreover, according to the definition of the acceleration factor in Equation (10), the acceleration factor can be maximized by making the drive torque τ i (t) of the individual gearbox have two states of ±τ max i . However, the drive torque τ i (t) obtained through such two state controls can cause chattering. In addition, the stepwise rapid change of drive torque can stimulate unmodeled dynamics such as the stiffness of the links and joints of the robot. This can lead to oscillatory robot motion, and the acceleration factor will also vary as the robot may not move as intended.
Therefore, in this study, the displacement of the robot drive component q(t), which is an optimization variable, is limited to have a sufficiently smooth profile that is continuous up to the second derivative. This ensures the resulting load torque τ i (t) of the robot drive component is continuous function of the time.
where k is the order of the B-spline, m is the number of control points, p j i is the j-th control point for interpolation of the displacement of the i-th joint, and B j,k (t) is the B-spline basis function.
According to the definition of B-spline curve, the knot sequence was defined as a sequence δ of m + k + 1 elements, and all elements in the sequence except for k + 1 overlapped elements at the beginning and end increase at equal intervals as the Equation (24) Appl. Sci. 2020, 10, 7459 When the displacement of the robot drive component is expressed using B-spline, the constraint Equations (21) and (22) on the displacement and speed of the optimization problem can be satisfied by setting four points of the B-spline control points as in Equations (25) and (26) [8].
Therefore, the optimal variable x is defined as in Equation (27) by augmenting m − 4 B-spline control points that describes the displacement of the drive components for each joint of the robot.
To summarize the above results, the constrained optimization problem of Equations (16)-(22) can be expressed as a finite-dimensional optimization problem of Equations (28)-(30) for the optimization variable x.
The linear inequality in Equation (29) represents the constraint in Equation (19) for the rotational speed of each joint in terms of the optimization variable x, and the inequality in Equation (30) represents the constraint condition in Equation (20) for the displacement of each joint.
The matrix A of Equation (29) is a block diagonal matrix and it is defined as a matrix having submatrix A i ∈ R (m−3)×(m−4) as elements. Submatrix A i is a band diagonal matrix and is defined as Equation (31).
The vector b in Equation (29) is an augmentation of vector b i ∈ R (m−3) , defined as in Equation (32). Then the deviation between adjacent control points of the B-spline was calculated by linear The vector c in the Equation (29) is an augmentation of vector c i ∈ R (m−3) , defined as in Equation (33), and the vector c i is derived by removing the first two and the last two elements from the knot sequence δ of Equation (24) and using the k -th difference of thus formed sequence where k is the B-spline order.
The vector x i , x i ∈ R (m−3) is defined by Equations (34) and (35) and includes parameters for the minimum and maximum angles of the individual joints.
The optimization tool utilized in this paper (python Scipy.optimize module) requires the objective function and the function for the constraint are entered separately, so if the constraint Equation (18) for the driving torque of each joint is reflected as a hard constraint, the dynamic equation must be evaluated in duplicate. In addition, it is necessary to implement a function with dimension as high as the number of sampling points to check whether the torque constraint is satisfied for each evaluation point. As a result, it may take longer computation time to solve the optimization problem. Therefore, we decide to reflect the constraint Equation (18) as a soft constraint by adding the quadratic penalty function of the vector ρ defined in Equation (36) to the objective function.
As a result, a constrained optimization problem with linear constraints on the control points of the B-Spline curve, which is an optimization variable, was established.

Problem Data
In order to analyze the dynamics of the 6-DoF industrial robot with the above algorithm, kinematic parameters such as the coordinate transformation matrices at home position M i = R i p i 0 1 and the screw vectors S i are necessary that appears in Equation (2). Here, R i represents a rotation matrix to convert a 3D spatial vector represented in the i-th link coordinates into the i − 1-th link coordinates. p i is a vector in which the origin of the i-th link coordinates is represented in the i − 1-th link coordinates. S i represents a screw vector to express the degree of freedom of the i-th joint. Table 1 shows the kinematic parameters for the industrial robot used in this study.  10, 7459 9 of 16 In the dynamics analysis with inward recursion, an inertial parameter expressed as

Index (i) Rotation (R i ) Translation (p i ) Screw Vector (S
is required. Here, I i is the rotational inertia at the center of gravity of the i-th link, r i is the position of the center of gravity of the i-th link expressed in the i-th link coordinates, and m i is the mass of the i -th link. Table 2 is the inertial parameters for the industrial robot used in this study. The inertial matrix J 7 for the object fixed to link 6 was converted into a representation of the coordinates of link 6 as  The limit on displacement, speed and torque of each joint is given in Table 3, where the rated torque τ rate for each joint is assumed to be same with the maximum torque τ max .  Table 4 presents information on the start and end postures of the robot joint displacement. The initial value of the B-spline control parameter, which is the optimization variable, was set to the motion profile moving at a constant velocity by applying the value obtained by dividing the start and end posture of the robot joint displacement by equal intervals.
The order of the B-spline describing the displacement of each joint was set to 3 to ensure C 2 continuity of the motion profile, and the number of spline control points for each joint was set to 10.

Implementation of the Methodology
The inverse dynamic algorithm presented in Equations (1)-(6) was implemented as a function using the kinematics and the dynamics parameter of the robot given in Tables 1 and 2, respectively. In order to evaluate the objective function given in Equation (28), it is necessary to calculate the acceleration factor in Equation (10) for the given motion profile. Therefore, the inverse dynamics algorithm must be called as much as the number of sampling points for numerical integration calculation in Equation (10). That means the inverse dynamics function must be executed t f ∆t × Number o f iterations times during an optimization, where ∆t is the time step (10 ms in this study) for the numerical integration.
As a result, quick processing of the inverse dynamics routine is critical for the computation time. To this end, we implemented the inverse dynamics routine with C++. The matrix operation required for the calculation of dynamics was implemented using the free software library Eigen (https://eigen.tuxfamily.org).
To solve the optimization problem formulated as in Equations (28)-(30), an optimization algorithm was implemented. The optimization algorithm initializes the matrix A and vectors b, c in Equation (29) using the problem data given in Tables 3 and 4. The bound x and x in Equation (30) are initialized using the problem data given in Table 3. The objective function for the optimization problem which is a part of the algorithm receives the B-spline control points as an input and evaluates the joint torques for each sampling points by invoking the inverse dynamic algorithm. Then, the objective function calculates the acceleration factor defined in Equation (10) for all joints of robot and returns AM-GM ratio of the acceleration factors.
The optimization algorithm was implemented in Python 3.5.2 environment. In particular, the interpolation module for SciPy library (version 1.2.1) was used to evaluate the B-spline value at the sampling points. Sequential Least Squares Programming (SLSQP) solver provided by the optimize package was also utilized to solve the constrained optimization problem (https://docs.scipy.org/doc/ scipy/reference/index.html).

Optimization Results
The optimization problem is solved for various cycle times on an Intel Core i9-8950HZ CPU 2.9 GHz personal computer with Linux (Ubuntu 16.04 LTS) operating system. The computation times to solve the optimization problems are given in Table 5.   Table 6 shows the result of optimization according to cycle time t f . The acceleration factors for the initial uniform motion and the optimal motion are compared for each cycle time. It can be observed that the initial acceleration factors with large deviations become uniform after the optimization. As the cycle time is decreased, the required speed to reach from the start position to the end position is increased, so the constraint on the speed is not satisfied if the cycle time is too short. In this study, optimization was performed by shortening the cycle time up to t f = 2 s.
The acceleration factors converge to 0.045 when the cycle time t f = 2, which means that about 22 (≈ 1/0.045) times longer test time is required for the system-level accelerated life test than that of component-level life test in which the drive components are tested under the rated conditions. Table 7 shows the mean, standard deviation, and coefficient of variation for the acceleration factor of each joint. Compared to the initial uniform motion, as a result of the optimization, it is confirmed that acceleration factor of each joint converges to the same value within a range of 0.1% or less based on the coefficient of variation.               Figures 5 and 6 represent the initial and the optimal motion profiles, respectively, for the case of t f = 2 s. Each joint begins to move at the initial angle and stops at the final angle, which are given in Table 4, at the end of the cycle time t f . In the initial motion, the speed is uniform except for the initial and final acceleration/deceleration of the cycle time, while in the optimal motion, the speed fluctuates significantly. In particular, for optimal motion, it can be seen that the speed of joint 1 and torque of joint 5 reach the lower limit in 0.6 s and remain within the bound in Table 3 for the rest of the time (Figure 6).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 16 Figures 5 and 6 represent the initial and the optimal motion profiles, respectively, for the case of = 2 s. Each joint begins to move at the initial angle and stops at the final angle, which are given in Table 4, at the end of the cycle time . In the initial motion, the speed is uniform except for the initial and final acceleration/deceleration of the cycle time, while in the optimal motion, the speed fluctuates significantly. In particular, for optimal motion, it can be seen that the speed of joint 1 and torque of joint 5 reach the lower limit in 0.6 s and remain within the bound in Table 3 for the rest of the time ( Figure 6). Since the spline order was set to 3, it was possible to obtain a continuous motion profile up to the second derivative of the displacement, and it was confirmed that the profile of the driving torque also has a continuous waveform.

Conclusions
In this study, a numerical method to design robot motions for system-level life test of drive components of robot is developed. The major results of the paper are summarized as follows.
In order to derive the robot motion that makes the acceleration factor of each joint equal, we formulated a robot motion optimization problem that minimizes the AM-GM ratio of the acceleration factor of each joint.
Starting from the initially linear motion, a global optimal robot motion was derived that make the objective function of the optimization problem converges to zero. The optimal robot motion was shown to be continuous and make the acceleration factor of each joint equal. Using the parameters of an actual 6-DoF industrial robot, global optimal robot motions were obtained for various cycle times and it was confirmed that the acceleration factor is inversely proportional to the cycle time. Since the spline order was set to 3, it was possible to obtain a continuous motion profile up to the second derivative of the displacement, and it was confirmed that the profile of the driving torque also has a continuous waveform.

Conclusions
In this study, a numerical method to design robot motions for system-level life test of drive components of robot is developed. The major results of the paper are summarized as follows.
In order to derive the robot motion that makes the acceleration factor of each joint equal, we formulated a robot motion optimization problem that minimizes the AM-GM ratio of the acceleration factor of each joint.
Starting from the initially linear motion, a global optimal robot motion was derived that make the objective function of the optimization problem converges to zero. The optimal robot motion was shown to be C 2 continuous and make the acceleration factor of each joint equal.
Using the parameters of an actual 6-DoF industrial robot, global optimal robot motions were obtained for various cycle times and it was confirmed that the acceleration factor is inversely proportional to the cycle time.
It is expected that the motion optimization method proposed in the present study can be utilized to design a test mode for performing a system-level accelerated life test in which the robot drive component composed of a servo motor and a gearbox are mounted on the robot.
The results obtained so far assume the cycle time and initial and final posture of the robot are the fixed data of the optimization problem and consider only the trajectory of joint variables as the optimization variables.
In future work, we will explore more efficient test modes by including the robot's cycle time, initial and final posture of robot, inertial parameter of the weight in the optimization variables. The effect of spline order and initial robot motion on the optimal solution will also be explored further in future studies.

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