Derivative Analysis Of Compensated Machine Coordinates Vectors

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

1. 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:

   1. Establish kinematics model

   2. Measurement

   3. Robots parameters identification

   4. Compensation

   Additionally, non-geometric 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 HA-06 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 70-90% 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.

2. 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 non-geometrical parameters in robot calibration procedure significantly improve the robot accuracy.

   Comparison between geometric and non-geometric 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 non-geometric 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 Paden-Kahan 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 3-5 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 G-code of machine coordinates should be produced in order to send it to the controller. That G-code 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 G-code 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 i-th machine

   coordinate in pose coordinates changes. Next, the differences

   machine configuration drivers and eventually violation of their limitations.

   F F x, i x

   (7)

3. 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 G-code 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 G-code, 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 G-code 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 G-code 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 G-code, 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 look-ahead 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 i-th 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 i-th

   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.

4. 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

1. 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

2. 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.131-146, ISBN: 3-86611-286-6, I-Tech, 2007

   zoom

   zoom

   Fig. 1 First, second and third derivative analysis of the machine coordinate 1 with respect to s

3. Nguyen, H-N., Zhou, J., Kang, H-J. A New Full Pose Measurement Method for Robot Calibration, Sensors 2013, 13, 9132-9147, DOI: 10.3390/s130709132, 2013

4. Santolaria, J., Conte, J., Ginéz, M. Laser tracker-based kinematic parameter calibration of industrial robots by improved CPA method and active retroreflector, The International Journal of Advanced Manufacturing Technology, 66(9-12), pp.2087-2106,

   DOI:10.1007/s00170-012-4484-6, Springer, 2013

5. Joubair, A., Bonev, I.A. Kinematic calibration of a six-axis serial robot using distance and sphere constraints, The International Journal of Advanced Manufacturing Technology, 77(1-4), pp.515-523,

   DOI:10.1007/s00170-014-6448-5, Springer, 2015

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

7. Zhou, J., Kang, H-J. A hybrid least-squares genetic algorithm-based algorithm for simultaneous identification of geometric and compliance errors in industrial robots, Advances in Mechanical Engineering 2015,

   Vol. 7(6) I-12, DOI:10.1177/1687814015590289, SAGE, 2015

8. Joubair, A., Nubiola, A., Bonev, I. Calibration Efficiency Analysis Based on Five Observability Indices and Two Calibration Models for a Six-Axis Industrial Robot, SAE International Journal of Aerospace 6, no. 2013-01-2117, pp. 161-168, DOI:10.4271/2013-01-2117, 2013

9. Kamali, K., Joubair, A. Bonev, I.A., Bigras, P. Elasto-geometrical Calibration of an Industrial Robot under Multidirectional External

   Loads Using a Laser Tracker, 2016 IEEE International Conference on Robotics and Automation (ICRA), pp. 4320-4327, DOI: 10.1109/ICRA.2016.7487630, IEEE, 2016

10. 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:I-II, DOI: 10.1177/1729881418755157, SAGE, 2018

11. Curry, H.B. The Method Of Steepest Descent For Non-linear Minimization problems, Quarterly of Applied Mathematics 2.3 (1944), pp. 258-261, 1944

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

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

14. Dimovski, I., Samak, S., Cvetkoska, D., Trompeska, M., Kochoski, F. Speed control in numeric controlled systems, Vol.39(LXV) No.1 2015 (49-63), ISSN 0351-336X, UDC: 519.713-74,

    , 2015

15. Sencer, B., Altintas, Y. & Croft, E. Feed optimization for five-axis CNC machine tools with drive constraints, International Journal of Machine Tools and Manufacture, 48. 733-745.

    DOI:10.1016/j.ijmachtools.2008.01.002,2008

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