 Open Access
 Authors : Samoil Samak , Vladimir Dukovski , Igor Dimovski
 Paper ID : IJERTV9IS110143
 Volume & Issue : Volume 09, Issue 11 (November 2020)
 Published (First Online): 19122020
 ISSN (Online) : 22780181
 Publisher Name : IJERT
 License: This work is licensed under a Creative Commons Attribution 4.0 International License
Derivative Analysis Of Compensated Machine Coordinates Vectors
Samoil Samak
Institute for Advanced Composites and Robotics IACR
Prilep, North Macedonia
Igor Dimovski
Institute for Advanced Composites and Robotics Prilep, North Macedonia
Vladimir Dukovski Faculty of Mechanical Engineering University St. Cyril and Methodius
Skopje, North Macedonia
AbstractBeside there are lot of papers about parametric robot calibration, there is lack of research focused on the last step of the parametric calibration procedure compensation. In this paper, an iterative numerical compensation algorithm is described in detail. In the case when the robots movement is controlled using its machine coordinates in order to compensate the robots errors, the issue is if such changes in machine coordinates eventually cause some unwanted changes or constraints violations in the feed or acceleration profile. Derivative analysis of compensated machine coordinates for robot motion typical in composite industry use case, showed there are not any violations of the constraints raised out of velocity, acceleration or jerk limitations, for all robots axes
KeywordsRobot calibration; compensation; derivative analysis

INTRODUCTION
In composite industry, the robot accuracy is one of the most important factors needed in the technological process. Any deviation from the given pose (position and orientation) of the robot tool could cause some product defect. In composite manufacturing, AFP (Automated fiber placement) or ATL (Automated tape laying) head, with very complex mechanical structure is used as a robot tool.
Most common defects in AFP/ATL technologies are gaps between fibers/tapes out of predefined tolerance, overlap of consecutive tapes, laminate wrinkling, local buckling of the first ply etc. More details about composite product defects, their sources and comprehensive analysis of such defects types are given in Jordaens & Steensels [1].
One can conclude, an improvement in the robot accuracy is one of the key factors in reducing composite product defects, despite it not the only factor. To avoid composite product defects, entire process should be controlled in real time, including defect detection and robot accuracy control.
Most common approach in robot accuracy improvement is parametric calibration. Most comprehensive review of robot accuracy improvement state of the art is given in Abderrahim et al [2]. They define the robot calibration as a process by which manipulator real parameters values affecting its
accuracy are established through direct or indirect measurement and used to improve its accuracy by modifying the positioning software.
Up to some modifications, the parametric calibration procedure follows these steps:

Establish kinematics model

Measurement

Robots parameters identification

Compensation
Additionally, nongeometric sources of errors as temperature, links stiffness, loads etc, in a similar manner could be included in the parametric calibration procedure.
Optotrac measurement system used to determine the tool position only and the Krypton K600 system with three linear CCD cameras developed by Metris company as well are described in [2]. Position accuracy is improved reducing the position error from 3.25mm to 0.29mm and the orientation error from 5.43mrad to 0.35mrad, for a robot with load of 158kg.
All the researches related to the robots parametric calibration mainly follow the 4 steps procedure described above. There are some variations in type of metrological equipment, choices of sample of measured poses in the robots workspace and whether nongeometric parameters are included or not.
Nguyen et al. [3] for example, have developed a model to calibrate 27 DH robot parameters, measuring set of positions that follow circle path. Their results are experimentally confirmed on Hyundai HA06 robot, reducing the average volumetric positional error from 3.657mm before the calibration, to 0.129mm after the calibration. Orientation error expressed by the Euler angles is reduced from [0.330o, 0.672o, 1.400o] in average, to [0.009o, 0.017o, 0.013o] after the calibration.
Similar research is conducted by Santolaria et al. [4], taking specific choice of measurement points made by Circle
Point Analysis approach. This approach advances are having clear geometric and physical meaning of the errors, since the robot tool is forced to describe circle movement, employing one particular robots axes at a time. They have used active target in order to achieve greater ranges measurements, but that increase the uncertainties. Some results can be found in
[4] for the experiments made on Kuka KR 5 sixx R650 robot. One can conclude beside most clear result in geometric representation of measured errors, improvement in total volumetric error reduction is significantly worse 33%40% compared to 7090% of improvement achieved by the models that perform parameters optimization in formal manner, without consideration of their physical meaning.Very accurate contact sensor in combination to specially designed artefact of 3 spheres is used to calibrate 25 robot parameters in the research presented by Joubair & Bonev [5]. The details of parametric calibration algorithm to identify 36 geometric parameters of KUKA 480 R3330 robot could be found in [6]. Only positional errors are taken into consideration and 2 levels of algorithm validation are provided using simulation and experimentally.
Zhou & Kang [7] go step beyond, including non geometrical parameters in the optimization procedure, using genetic algorithm. They emphasize the stiffness of the robots
orientation accuracy are obtained, reducing the total orientation error from 0.147o to 0.095o in average.
Making literature review, one can notice less information are provided for the step for in typical parametric calibration procedure compensation. Usually, calibrated parameters are sent to the robot controller and there is lack of technical details how are they used to compensate the robots errors. This issue is especially important from robots kinematic and dynamic analysis point of view, in the case when the robots movement is controlled using its machine coordinates in order to compensate the robots errors. Interpolation also has some impact on the robot accuracy and some fluctuations could eventually lead to unwanted jumps in the velocity or acceleration profiles, for some of the robot axes.
In the following section, our iterative numerical compensation algorithm is elaborated in detail. In the section 3, derivative analysis of the machine coordinates is provided. Conclusions are given in the last section.


ITERATIVE NUMERICAL COMPENSATION ALGORITHM
Let Pdesired be a given pose vector of the 6R robot tool and its scalar components with respect to the reference coordinate system (RCS) are:
RCS
RCS
links as key reason why geometric parameters calibration does not provide maximal accuracy improvement across entire workspace, especially in a case of complex tool as significant
P
desired
X , Y , Z , A, B, CT
(1)
load. Experimental verification shows average positional error reduction to 0.091mm. Compared to results in [3], including nongeometrical parameters in robot calibration procedure significantly improve the robot accuracy.
Comparison between geometric and nongeometric robot calibration is given in Joubair er al. [8]. The experimental data
Te first three coordinates of the desired pose Pdesired refer
to the position of tool center point (TCP) in millimeters, and the last three coordinates are Euler angles in degrees, and they refer to the tool orientation with respect to RCS. This format of desired pose in transformed to homogeneous 4×4 matrix for the purpose of easy calculations:
given there refer to Fanuc LR Mate 200 iC robot with both of
nx lx ox
px
approaches and using more than 7000 measurement points. Average positional error before calibration is 0.622mm. After the geometric calibration it is reduced to 0.250mm and after including stiffness of the last 5 robots axes, it is reduced to
n l o
T y y y
n l o
z z z
py
p
z
(2)
0.142mm. Tool orientation is not taken into consideration and an important note is given that only 200 measurement points are sufficient to identify all the robot parameters. More details how to include nongeometric parameters in the mathematical model are given by Kamil et al. [9]. In this research, different loads are used in order to experiment the impact of stiffness parameters. Experimental data refer to ABB IRB 1600 robot and the maximal positional error is 0.960mm if stiffness parameters are used, compared to 2.571mm if just geometric parameters are used.
In our research, we have designed parametric robot
0 0 0 1
In order to determine appropriate machine coordinates for the given pose, the inverse kinematics algorithm based on the PadenKahan subproblems could be applied. Details of algorithmic approach of the generalized form of these subproblems are given in Dimovski et al [10]. Using this approach, 8 solutions of machine coordinates are usually determined, so trajectory planning algorithm is used to choose the most optimal one. For 6R industrial robot used in this
T
T
research, 6 machine coordinates respective to 6 rotational axes are obtained:
calibration algorithm including 30 geometric parameters based on the kinematic model built as screw theory model [10] and 5 stiffness parameters. This approach is experimentally
ideal
1, i
, 2, i
, 3, i
, 4, i
, 5, i
, 6, i
(3)
confirmed on machine configuration consists of KUKA KR500 robot with ATL head manufactured by Mikrosam company. The average positional error is reduced from 2.040mm before the calibration, to 0.580mm after the calibration. Orientation deviations are included in the optimization procedure as well. Improvements in the
There are at least two reasons why this approach is idealized and mostly theoretical. First, nominal robot parameters are user, and second assumptions like last 3 axes of rotation are concurrent they intersect at a single point in 3D space. In reality, there are always some deviations of such idealized assumptions. The point is such deviations are
quantified in the calibration procedure, so calibrated parameters should be used instead of nominal ones.
2
compensated
1
compensated
J 1 F
x
(8)
No matter what kind of calibration algorithm is applied, the new set of calibrated parameters is obtained, and they
The result of the first iteration
2
compensated
is used as
differ from the nominal parameters. Also, these calibrated values of the parameters cannot be simply substituted in the same inverse kinematics procedure, since changing the parameters values, the assumptions of parallelism, orthogonally and concurrency are not valid anymore.
Because of that, iterative numerical compensation
input in the second iteration. Completely the same procedure is applied until the objective function defined as a square of the magnitude of the vector Fx become small enough, within the tolerance . That means iteration procedure stops when the condition:
algorithm is designed in order to calculate compensated machine coordinates appropriate to the given desired pose
2
F
x
(9)
Pdesired, using the new set of calibrated values of the robot parameters.
, , , , , T (4)
is satisfied.
In this research, described iterative algorithm is implemented using the tolerance value of 1015 and
compensated 1, c 2, c 3, c 4, c 5, c 6, c
iterative step value of 108 . In most cases the iterative
This algorithm is based on the Newton method. This approachs fundamentals are given by Curry [11]. In the contemporary applied mathematics literature, more specifically papers refer to optimization problems, some details of the Newton method, its application and variety of modifications and implementations could be found, like in [12], [13].
The desired pose coordinates Pdesired are input in this algorithm, as well as the ideal machine coordinates ideal, used as initial guess. The algorithm actually iteratively makes deviation reduction in order to get as closer as possible to the desired pose Pdesired, in the pose space, to some predefined tolerance that also comes as input in the algorithm. One additional input value is needed iterative step as a measure of machine coordinates changes in each iteration.
In the first iteration, machine coordinates are set to the initial guess values:
numerical compensation algorithm converges after 35 iterations.
Method for calculating the compensated machine coordinates after parametric robot calibration procedure, for given particular point Pdesired in the pose space is described in detail above. In practice, when 6R industrial robot, or even more complex machine configuration is used in composite tape laying is controlled by machine coordinates, usually an array of desired poses comes as input and Gcode of machine coordinates should be produced in order to send it to the controller. That Gcode is produced in some offline procedure such that for every desired pose, the described iterative numerical compensation algorithm is called in order to calculate appropriate compensated machine coordinates
compensated, based on calibrated values of the robot parameters, no matter what kind of parametric calibration procedure is performed to derive such calibrated robot parameters.
Beside the given desired poses, additionally machine
1
compensated
ideal
(5)
configuration trajectory planning is very important in composite manufacture, since large number of machine
T
T
The forward kinematics algorithm with calibrated values of the parameters is called in the next step and its output is used to determine deviation from the desired pose Pdesired. Such deviations are calculated and transformed into the pose space in order to obtain 6 scalar components of the deviations vector:
coordinate vectors need to be interpolated through the Gcode points and tool path has to be controlled and highly accurate along all the path. In online interpolation procedure, but in offline as well, the iterative numerical compensation algorithm calculations need to be performed very fast, and in the online case in real time.
F
x X
, ,
Y Z
, A, B
, C
(6)
Challenge is additional check if such changes in the
machine coordinates due to compensation cause unwanted changes in feed, acceleration or jerk for each of the robot or
The vector Fx,i determines the contribution of ith machine
coordinate in pose coordinates changes. Next, the differences
machine configuration drivers and eventually violation of their limitations.
F F x, i x
(7)

DERIVATIVE ANALYSIS OF THE MACHINE
are calculated for every i = 1,2,3,4,5,6. These 6 vectors with length 6 are organized as 6×6 matrix and that matrix divided by the iterative step is actually Jacobian matrix J.
COORDINATES AS A FUNTION OF TOTAL AXES DISPLACEMENT
In order to check if some unwanted changes and violation of feed, acceleration or jerk constraints occurs, the array of
The machine coordinates
1
compensated
are then
compensated coordinates vectors is analyzed against appropriate uncompensated coordinates vectors.
corrected by Jacobian inverse multiplied by the vector Fx:
Tool path defined as line in 3D physical space, with length of approximately 1500mm is the use case analyzed in this research. Since the tool path is straight line, theoretically it could be defined using Gcode of only 2 points. But, in this case, the interpolation between these points is completely left to the controller and the tool orientation is not controlled as well. That is inadmissible, especially in the case of robot usage in composite tapes laying, since then following the predefined position and orientation of the tool along all the path is extremely important. In the practice, additional number of poses is offline programmed, usually taking uniform Euclidian distance travel of the tool center point.
In this use case, the straight line is divided to 297 points at a uniform distance of 5mm and positions and orientations are determined for each of these points as in equation (1). Taking the nominal robot parameters and calling the standard Paden
function with respect to the time t. Even more, feed, acceleration and jerk constraints for each machine configuration axis is expressed with respect to total vector displacement s, allowing trajectory planning, optimizing the feed and in this case, allowing a simple analysis if there are any constraints violations, before trajectory planning is done.
More details how such constraints could be derived are given by Altintas & Kaan [16].
Starting from 297 vectors of machine coordinates given with appropriate Gcode, total vector displacement of the robot axes is calculated using Euclidian metric in 6D machine space. The total axes vector displacement s is uniformly divided to 2900 parts, and for every appropriate value of s, interpolated values of each machine coordinate is obtained by linear interpolation:
Kahan inverse kinematics procedure,
vector of
qs s, s, s, s, s, sT
(10)
ideal
1 2 3 4 5 6
machine coordinates is calculated, for each of 297 poses. This
The velocity of the axis i is denoted by (i=1,2,3,4,5,6):
way, the Gcode of uncompensated machine coordinates is
di
(11)
created.
The same 297 poses are used calling the iterative
i dt
numerical compensation algorithm, described in the previous section, in order to obtain another Gcode of compensated machine coordinates, concerned to the same path, the tool should follow.
Tool poses refer to the Pose Space and Euclidian distances used to generate such uniformly distributed poses along the tool path refer to the physical 3D space where the robot motion is done. That does not mean each of the robot axes
The acceleration of the axis i is denoted by:
i
i
d 2
i dt 2
The jerk of the axis i is denoted by:
i
i
d 3
i dt 3
(12)
(13)
would perform uniform displacement. Controlling the robot motion sending the machine coordinates to the robot controller, the focus of interest is moved to the Machine Space. In the case of 6R industrial robot used in this research it is 6D machine space. Beside between any 2 successive points in the Gcode, the tool always should travel 5mm displacement, in general, every robot axis should perform different displacement if any. As a consequence, velocities of the robot axes are different in different point of time. In this case, trajectory planning is left to the controller and it applies its own algorithms to interpolate several values taking into account to not violate the machine configuration dynamics.
In order to take full control of each robot axis motion, some algorithm of the lookahead class of algorithms should be applied, so after the complete trajectory planning every machine coordinate would be expressed as function of time t.
Analysis of the performances for some of these algorithms can
Consequently, the vectors of velocity, acceleration and jerk are denoted by q s , q s and q s respectively.
Evaluation of these functions can be done only after the trajectory planning. In fact, they are used in trajectory planning optimization – machine configuration motion program, minimizing the total time needed to pass entire path and keeping all feed, acceleration and jerk constraints satisfied.
Equation (10) determines an array of values for ith machine coordinate (i=1,2,3,4,5,6), for each of 2900 values of total vector displacement s of the robot axes.
Numerical procedure for approximation of the first, second and third derivative with respect to the total vector displacement s are applied. These derivatives are denoted as follows:
be found in [14]. Application of such trajectory planning
di
(14)
algorithm to express the robot machine coordinates with
i,s ds
respect to time t is out of the scope of this particular research.
One of the most important researches in the area of machine configuration trajectory planning, considering the
2
d
d
, i
i,ss ds2
d
d
3
i
(15)
feed optimization is elaborated by Sencer et al. [15]. Nomenclature used in this paper is especially useful and it
i,sss
ds3
(16)
allows making clear distinction between trajectory planning as a function with respect to total vector displacement s of the axes in the machine space and optimization of the velocities of each of the axes expressing the axis displacement as a
This way, vectors containing the first, second and third derivative of the 6 robot axes positions with respect to total vector displacement s are obtained:
q s
s,
s,…,
sT
(17)
each robot axis, as a function with respect to total vector
s
qss
1,s
1,ss
1,ss
s
2,s
s,2,ss
6,s
s,…,6,ss
sT
displacement s are determined by:
(18)
T
T
qsss s 1,sss s,2,sss s,…,6,sss s
(19)
sA , maxs
a
(27)
max i 1,2,3,4,5,6
max i 1,2,3,4,5,6
qi,ss s
qi,ss s
1
1
1
1
If the maximal values of each axis velocity are given, the
normalized partial derivatives could be expressed:
sJ , maxs 3
j
(28)
qv
i,s
s
(20)
max i 1,2,3,4,5,6
qi,sss s
i,s
vi,max
Visual analysis of the first, second and third derivative of
The velocity (feed) maximal allowed values, for each robot axis, for particular value of total vector displacement s,
are determined by:
the machine coordinate 1, referred to A1 robots axis is shown on Fig. 1. First four diagrams show the displacement 1 s and its derivatives 1,s s , 1,ss s and 1,sss s with respect to
1 s, taking
s 0.4;1.5rad . Taking wider range of the total
i,s
i,s
sV maxs
max
i 1,2,3,4,5,6
qv s
(21)
vector displacement s allows an analysis of the changes of appropriate functions, how frequent such changes happen, are
Derivation of the equation (21), as well as acceleration and jerk constraints derivation details could be found in [6].
The equation (20) actually determines the first derivatives constraints of the functions i,s s , raised out of velocity
limitations. The first derivatives have to satisfy the constraints as well raised out of the acceleration and jerk limittions.
s
there any unwanted fluctuations and if there is some significant difference between compensated and uncompensated values. If there are some violations of the limitations raised out of velocity, acceleration and jerk of the robots first axis, could be noticed as well.
Blue graphs refer to the uncompensated values and the red ones to the compensated values of the machine coordinate 1 of the robots axis A1, for all diagrams.
q a
q a
a i,s i,s
i,max
(22)
One can conclude there are not significant deviations between uncompensated and compensated values; especially
q
q
j i,s
i,s s
ji,max
(23)
there are not such deviations that may cause limitations violation. The red graphs are slightly shifted compared to the blue ones, without significant increase or decrease of the
In total, there are 3 constraints for the first derivatives of the functions i,s s . In the equations (22) and (23), notations
maximal or minimal values in the observed interval.
First diagram simply shows the machine coordinate 1 s
ai,max
and
ji,max
are acceleration and jerk limitations for ith
as a function with respect to s.
axis respectively.
Similarly, the second derivative of the functions s
Second diagram shows the first derivative 1,s s and as
i,s
with respect to the total vector displacement of the robot axes has to satisfy the following 2 constraints raised out of acceleration and jerk limitations respectively, for each robot axis:
well the constraints raised out of velocity limitations of the appropriate robots axis, shown as dotted green line; the constraints raised out of acceleration limitations of the appropriate robots axis, shown as dotted magenta line; and the constraints raised out of jerk limitations of the appropriate
qa
i,ss
s
robots axis, shown as dotted cyan line. One can conclude there are not violations of all the constraints at all.
i,ss
ai,max
Third diagram shows the second derivative 1,ss s and as
q
q
j
i,ss
i,ss s
j
(24)
well the constraints raised out of acceleration and jerk limitations of the robots axis A1, shown in similar manner
i,max
(25)
using magenta and cyan dotted line respectively.
Finally, the fourth diagram shows the third derivative
The third derivative of the functions
i,s
s with respect to
1,sss s and the constraints raised out of jerk limitations of the
the total vector displacement of the robot axes has to satisfy 1 constraint raised out of jerk limitations, for each robot axis:
robots axis A1, shown as dotted cyan line.
This analysis shows there are not limitations violations
q
q
j
i,sss
i,sss
j
s
(26)
neither for first, second or third derivative on the observed interval.
i,max
According to that, the maximal allowed feed values, satisfying the maximal acceleration and jerk limitations, for
In order to capture the more frequent changes of the presented function, the interval s 0.710;0.725rad is used
in last four diagrams to zoom the diagrams from the previous four diagrams.

CONCLUSION
Parametric calibration algorithm is designed and is experimentally confirmed on machine configuration consists of KUKA KR500 robot with ATL head manufactured by Mikrosam company. As result, the average positional error is reduced from 2.040mm before the calibration, to 0.580mm after the calibration and the total orientation error is reduced as well from 0.147o to 0.095o in average.
In order to include identified robot parameters in the robot control, using machine coordinates, iterative numerical compensation algorithm is designed and explained in detail.
The issue if some unwanted changes would appear as result of machine coordinates changes, derivative analysis of the machine coordinates as a function of total axes displacement is provided. As use case, typical robot motion in composite industry is used and derivative analysis showed there are not any constraints violations neither related to velocity, acceleration or jerk limitations of the robot axes.
REFERENCES

Jordaens, A., Steensels, T. Formation of defects in flat laminates during automatic tape laying (framework of a masters thesis), Faculty of engineering technology, Leuven, Belgium, 2015

Abderrahim, M., Khamis, A., Garrido, S., Moreno, L. Accuracy and Calibration Issues of Industrial Manipulators chapter in Low, K.H. (ed) Industrial Robotics – Programming, Simulation and Applications, pp.131146, ISBN: 3866112866, ITech, 2007
zoom
zoom
Fig. 1 First, second and third derivative analysis of the machine coordinate 1 with respect to s

Nguyen, HN., Zhou, J., Kang, HJ. A New Full Pose Measurement Method for Robot Calibration, Sensors 2013, 13, 91329147, DOI: 10.3390/s130709132, 2013

Santolaria, J., Conte, J., GinÃ©z, M. Laser trackerbased kinematic parameter calibration of industrial robots by improved CPA method and active retroreflector, The International Journal of Advanced Manufacturing Technology, 66(912), pp.20872106,
DOI:10.1007/s0017001244846, Springer, 2013

Joubair, A., Bonev, I.A. Kinematic calibration of a sixaxis serial robot using distance and sphere constraints, The International Journal of Advanced Manufacturing Technology, 77(14), pp.515523,
DOI:10.1007/s0017001464485, Springer, 2015

Samak, S., Dimovski, I., Trompeska, M., Hristoski, M., Kochoski, F., Dukovski, V. Computerbased simulation and validation of robot accuracy improvement method and its verification in robot calibration procedure, ETAI 2018 International Conference, ETAI Society of Macedonia, 2018

Zhou, J., Kang, HJ. A hybrid leastsquares genetic algorithmbased algorithm for simultaneous identification of geometric and compliance errors in industrial robots, Advances in Mechanical Engineering 2015,
Vol. 7(6) I12, DOI:10.1177/1687814015590289, SAGE, 2015

Joubair, A., Nubiola, A., Bonev, I. Calibration Efficiency Analysis Based on Five Observability Indices and Two Calibration Models for a SixAxis Industrial Robot, SAE International Journal of Aerospace 6, no. 2013012117, pp. 161168, DOI:10.4271/2013012117, 2013

Kamali, K., Joubair, A. Bonev, I.A., Bigras, P. Elastogeometrical Calibration of an Industrial Robot under Multidirectional External
Loads Using a Laser Tracker, 2016 IEEE International Conference on Robotics and Automation (ICRA), pp. 43204327, DOI: 10.1109/ICRA.2016.7487630, IEEE, 2016

Dimovski, I., Trompeska, M., Samak, S., Dukovski, V., Cvetkoska, D. Algorithmic approach to geometric solution of generalized Paden Kahan subproblem and its extension, International Journal of Advanced Robotic Systems, 2018:III, DOI: 10.1177/1729881418755157, SAGE, 2018

Curry, H.B. The Method Of Steepest Descent For Nonlinear Minimization problems, Quarterly of Applied Mathematics 2.3 (1944), pp. 258261, 1944

Nestorov, Yu. Introductory Lectures on Convex Optimization: A Basic Course, Springer, 2004

LemarÃ©chal, C. Cauchy and the gradient method", Doc Math Extra, 251 (2012): 254, 2012

Dimovski, I., Samak, S., Cvetkoska, D., Trompeska, M., Kochoski, F. Speed control in numeric controlled systems, Vol.39(LXV) No.1 2015 (4963), ISSN 0351336X, UDC: 519.71374,
, 2015

Sencer, B., Altintas, Y. & Croft, E. Feed optimization for fiveaxis CNC machine tools with drive constraints, International Journal of Machine Tools and Manufacture, 48. 733745.
DOI:10.1016/j.ijmachtools.2008.01.002,2008

Altintas, Y., and Kaan, E. "Feedrate optimization for spline interpolation in high speed machine tools", CIRP Annals 52.1 (2003): 297302, 2003