Robot dynamic model: freudenstein-based optimal trajectory and parameter identification

Abstract A thorough understanding of an industrial robot’s dynamic model is critical for practical robotic applications. An effective dynamic model is required for optimal controller design and trajectory planning. Robot manufacturers only provide kinematic data, which can only guarantee a certain level of positioner accuracy. The design of the trajectory-planning scheme, on the other hand, necessitates a thorough understanding of its dynamic features. The identification of dynamic parameters involves several procedures. This study used the Euler-Lagrangian equation to derive the robot dynamic model in its canonical form. To excite each link of the irb1600 robot industrial robot while avoiding displacement, velocity, and acceleration discontinuities at the start and endpoints, a Freudenstein 1-3-5 trajectory based on Fourier series expansion was used. The dynamic parameters were determined using the nonlinear least-squares approach based on the Levenberg-Marquart equation. The Savitzky-Golay smoothing filters improved the identification method by decreasing system noise.


ABOUT THE AUTHOR
The authors are members of the Gibela Research Chair in Manufacturing and Skills Development at Tshwane University of Technology, South Africa. The research chair is led by Prof. K Mpofu. The team's research cuts across the following niche areas: Expert Systems, Artificial Intelligence, Robotics in Manufacturing and Mechatronics, Virtual Commissioning Systems, and ROS-I applications are all being used to improve manufacturing. The team is involved in numerous projects: merSETA Localisation Project For Incubation And Skills Development; Virtual Commissioning in Transport Manufacturing Enterprises; Demanufacturing Research focusing on the locomotive; Smart Mobility Systems; and the Smart Factory Project.

PUBLIC INTEREST STATEMENT
A thorough understanding of an industrial robot's dynamic model is critical for practical robotic applications. An effective dynamic model is required for optimal controller design and trajectory planning. This study used the Euler-Lagrangian equation to derive the robot dynamic model in its canonical form. To excite each link of the irb1600 robot industrial robot while avoiding displacement, velocity, and acceleration discontinuities at the start and endpoints, a Freudenstein 1-3-5 trajectory based on Fourier series expansion was used. The dynamic parameters were determined using the nonlinear least-squares approach based on the Levenberg-Marquart equation. Therefore, a deep understanding and the identification of the base parameters can be used to develop an advanced path planning scheme and model-based controller. The requirement to identify the dynamic robot model stems from a lack of understanding of the robot model characteristics, which results in an inaccurate dynamic model and realistic simulation with practical implications.

Introduction
Robotic applications' usefulness in production is contingent on the accuracy of its tip, owing to the need for precise tracking at high speeds. Numerous variables influence the accuracy and dynamic performance of industrial robots: joint friction, tracking errors of the controller, tolerance in the machining, drive nonlinearities and flexibilities in the drives and joints (Carvalho Bittencourt, 2014;Gautier & Briot, 2013;Hardeman, 2008;Madsen, 2020;Pagani et al., 2020;Sweet & Good, 1985;Villagrossi et al., 2014). Therefore, a deep understanding and the identification of the base parameters can be used to develop an advanced path planning scheme and model-based controller. In robot manipulators, identification is divided into two categories: identifying the robot's kinematic structure and detecting the inertial characteristics of the linkages, actuators, and unknown load the robot carries (Kozlowski, 2012). To effectively deploy robotic techniques in manufacturing, as outlined in (Ogbemhe & Mpofu, 2015;Ogbemhe et al., 2019), a good understanding of the underlying kinematic and dynamic characteristics is required. System identification is the general field concerned with identifying models from measurements (Siciliano & Khatib, 2008). There are two models for system identification: parametric and non-parametric models. Non-parametric models are for complex systems where a few lumped parameters represent the system, whereas parametric models describe many parameters. Parametric models are ideal for manipulators since they are controllable and feature human-made components. The requirement to identify the dynamic robot model stems from a lack of understanding of the robot model characteristics, which results in an inaccurate dynamic model and realistic simulation with practical implications.
There are various works of literature on the estimate of industrial robot dynamic parameters (Ding et al., 2015;Hao et al., 2021;Liu et al., 2021;Madsen, 2020;Pagani et al., 2020); these methods are classed as offline and online estimation strategies (see, Figure 1).
Physical experimentation, computer-aided design (CAD), and parameter identification are examples of offline estimation methods. The online estimation technique uses advanced adaptive control and neural network algorithms. On the other hand, the online estimating technique employs sophisticated adaptive control and neural network technologies. For example, Shang and Cong (2015) demonstrated that parallel robots could enhance tracking accuracy by identifying and calibrating two degrees of freedom. These parameters may now be calculated from robot CAD models by leveraging the capabilities of contemporary CAD software, which is continually expanding. Several researchers have expressed an interest in identifying dynamic parameters. The authors of Gao et al. (2018) employed singular value decomposition (SVD) to identify linearly dependent robotic structural parameters and created an identification strategy that avoids being trapped in a local minimum. Gautier and Briot (2013) studied the use of the Power Identification Model to detect robot inertial parameters by constructing a block-triangular observation matrix using optimised trajectories. Using a non-gradient technique in dynamic parameter identification is becoming more prevalent.
Developing an ideal robot excitation trajectory requires nonlinear optimisation with motion constraints: limitations on the joint angle, velocity, acceleration, and the robot end-effector pose in Cartesian space to avoid colliding with the environment. Numerous ways have been advanced; they all utilised a different method for parameterising the trajectory. For example, in Ding et al. (2015), the authors used an artificial bee colony technique to discover unknown parameters. They validated their strategy using torque prediction accuracy. The method of excitation utilised to excite the robot trajectory is critical in establishing its dynamic properties. Polynomial Huang et al. However, in some cases, these approaches cannot minimise the amplitude of the acceleration profile to reduce effort on the load caused by the vibrationtional effects of the mechanical structure. Also, a discontinuous acceleration profile lowers the peak of acceleration but may induce oscillations and vibrations due to the corresponding discontinuity of the inertial forces. Furthermore, using Fourier series expansion trajectories can lower the amplitude of the acceleration profile while ensuring the continuity of the position profile and its derivatives up to a specific order (Biagiotti & Melchiorri, 2008). A robust identification model requires optimality criteria during excitation. Kröger (2010), Jin and Gans (2015) proposed Hadamard's inequality to optimise excitation trajectory using an understandable and straightforward criterion. Villagrossi et al. (2014) offer an index for analysing correlation influence among critical parameters under kinematic coupling limitations to determine optimal excitation patterns. Ding et al. (2015) used a fifth-order Fourier series as an excitation trajectory by minimising the minimal condition number. In a related effort, an ideal excitation trajectory for a DLR Light-Weight Robot was described as a restricted optimisation problem using parameterised B-splines (Rackl et al., 2012). The weighted square (Calanca et al., 2010) and least square methods (Rackl et al., 2012) have been widely employed for parameter estimation in industrial robotics. However, while the least square technique produces a more straightforward calculation, it is sensitive to measurement noise and necessitates data pre-processing. The weighted least squares approach, on the other hand, handles issues that are sensitive to noise.
This study derives the robot dynamic model in its canonical form. The industrial robot is excited using a Freudenstein 1-3-5 based optimal trajectory and parameter identification. The use of more Fourier series expansion terms in the Freudenstein 1-3-5 results in profiles with slower acceleration but higher frequency components during excitation. The remainder of the paper is organised as follows. Section 2 discusses the formulation of the dynamic robot model. Sections 3 and 4 deal with the dynamic identification model and experimental validation.

Robot dynamic model
The robotic manipulator is defined by each link's kinematic and dynamic parameters. The dynamic model is derived using the theory of generalised coordinates and a scalar function known as the Lagrangian; it contains the difference between the mechanical system's kinetic and potential energy. As a result, an accurate dynamic model is required for robot identification. Figure 2a illustrates the ABB irb1600 schematic and coordinates frame. Figure 2b displays the structure of an industrial robot. The Denavit Hartenberg (DH) specifications of the robot are shown in Table 1.
This study aims to present a parameter identification approach, which can correct for the uncertainties in the robotic system during excitation. First, it is mentioned that the examined robot is with 6 degrees of freedom (DOF), which is coupled in series. Moreover, the primary objective of this work is to present a Freudenstein 1-3-5 trajectory, built on the foundation of Fourier series expansion, which was used to excite each link of the industrial robot. A slower acceleration rate is achieved in the Freudenstein 1-3-5 profiles due to the increased usage of Fourier series expansion terms.
The dynamic equation of the industrial robot can be derived from the Lagrange-Euler formulation. It is written in matrix form as  τ f t ð Þ 2 R n is the friction forces and τ t ð Þ 2 R n represents the joint torque vector that is the input to the system. The friction force is modelled as where the f v and f c are constants n � n diagonal matrices representing viscous and Coulomb friction parameters, respectively, and sgn � ð Þ is the sign function. The dynamics of the entire sixaxis robot can be modelled by: The dynamic equation for the second link is derived in terms of position and time derivatives in joint space. It takes the following general form in (5) with inertia parameters described in Table A1.

Forces acting on the third joint
The dynamic equation for the third link is derived in terms of position and time derivatives in joint space. It takes the following general form in (6) with inertia parameters described in Table A1.

Forces acting on the fourth joint
The dynamic equation for the fourth link is derived in terms of position and time derivatives in joint space. It takes the following general form in (7) with inertia parameters described in Table A1.

Forces acting on the fifth joint
The dynamic equation for the fifth link is derived in terms of position and time derivatives in joint space. It takes the following general form in (8) with inertia parameters described in Table A1.

Dynamic identification model
The complexity of the robot's dynamic equations makes estimating its parameters problematic. On the other hand, the equations of motion are linear in terms of the inertia parameters.
(1) describes the mathematical model in joint space, which may be rewritten as a linearly parameterised form using the DH convention as indicated in (10): where Y θ; _ θ; € θ À � 2 R n�Ns is called the regressor matrix and X s 2 R Ns�1 is the standard parameter vector. In summary, Figure 3 illustrates the parameter identification method's operations; meanwhile, the suggested approximation identification model does not affect the inertia and joint viscous parameters in X dependency on them. Additionally, analogous to the rigid body case, the exciting trajectories may be optimised across the robot's workspace to mitigate the influence of small perturbations and ensure identification accuracy.
The six components of the inertia matrix, link masses viscous, and column friction coefficients are among the 13 standard parameters for a rigid industrial robot through each link and joint (Jin & Gans, 2015). Furthermore, an excitation reference trajectory is continually used to stimulate the industrial robot. Therefore, we use a Freudenstein 1-3-5 based optimal periodic trajectory in this study. The use of the Freudenstein 1-3-5 based trajectory to excite the industrial robot is motivated by the following: • It may be advantageous to reduce the amplitude of the acceleration profile to prevent load efforts caused by inertial forces or mechanical structural vibrations.
• A discontinuous acceleration profile lowers the peak of acceleration but may induce vibrations due to the corresponding discontinuity of the inertial forces.
• On the other hand, the cycloidal trajectory has a low harmonic content but has a greater acceleration value.

Nonlinear least square
The nonlinear least square (NLS) problem described in this manuscript is modelled as (Gavin, 2019): where ŷ t i ; p ð Þ is a model function of an independent variable t and a vector of n parameters p to a set of k data points t i ; y i ð Þ. The variable W is the weighting matrix and σ yi is the measurement error for the datum y t i ð Þ.
Because it does not need as many equations as the number of variables, the Levenberg-Marquart technique is employed to construct our nonlinear least square algorithm rather than the Trust-region-reflection approach. Furthermore, the Levenberg-Marquart uses two numerical minimisation techniques: gradient descent and Gauss-Newton, by adaptively altering parameter updates between the two schemes using the update mechanism in the Levenberg-Marquart (12), where λ is the damping parameter and J is used to denote @ŷ @p . @ŷ=@p

Trajectory parametrisation and optimisation
A Freudenstein 1-3-5 trajectory is used to determine the excitation of each link; it is based on Fourier series expansion. The Freudenstein 1-3-5 trajectory was chosen because it minimises the amplitude of the acceleration profile, hence avoiding mechanical structural vibrations. This trajectory is expressed as (Biagiotti & Melchiorri, 2008): where h ¼ θ f À θ 0 ; θ 0 is the initial displacement; θ f is the final displacement; T is the duration of the trajectory; t 0 is the initial time instant; t 1 is the final time instant; α is 0.9438.
The problem of finding an excitation trajectory θ � t ð Þ is formulated as where J is an objective function to be determined afterwards, and θ max i°ð Þ, _ θ max i°s À 1 À � , and € θ max i°s À 2 À � denote bounds on joint position, velocity and acceleration, respectively. The constraints in (17) are necessary to avoid discontinuities in displacement, velocity and acceleration at the start and endpoints in case θ i t 1 ð Þ�0, _ θ i t 1 ð Þ�0, and € θ i t 1 ð Þ�0.

Testbed and signal processing
The ABB irb1600 industrial manipulator has been employed in this simulation experiment depicted in Figure 2(a), and its associated link frames are presented in Figure 2(b) (b). The connection parameters of the robot are shown in Table 1. The ABB irb1600 has a serial structure with six rotation joints, as shown in Figure 3, making it ideal for glueing, arc welding, machine tending, material handling, deburring, and grinding tasks. In addition, the robot comes with the IRC5 controller, control software, and FlexPendant. However, the data gathered during identification, such as displacement, velocity, and acceleration, is prone to noise. Therefore, improving the accuracy of the parameter identification results requires reducing noise on (11) and (12). This study uses Savitzky-Golay smoothing filters to eliminate a noisy signal with a wide frequency range. The Savitzky-Golay filters (SGF) outperform traditional Finite Impulse Response (FIR) filters (Schafer, 2011), which are prone to filtering high-frequency information with the disturbance. Thus, they can help reduce the leastsquares error when fitting a polynomial to noisy data frames. The sgolayfit function in MATLAB 2020a is used to implement the Savitzky-Golay filters (Schafer, 2011). The PC's specs include an Intel (R) Core(TM) i3-3110 M processor, 12GB of RAM, and a 1TB hard drive. Figure 4 shows that the reference trajectory for all joints is smooth, reducing the robot arm's shock and impact and ensuring the smooth operation of the ABB irb1600 robot. The black line represents the displacement after nonlinear least-squares identification with a disturbance. In contrast, the red line represents the displacement after applying the Savitzky-Golay smoothing filter. Furthermore, there is a significant disturbance in the third joint during excitations compared to the other joints. Figure 5 demonstrates that the velocity for all joints is smooth, minimising the robot arm's shock and impact and assuring the smooth functioning of the ABB irb1600 robot. The black line shows the velocity following nonlinear least-squares identification with a disturbance. In contrast, the blue line indicates the velocity following the Savitzky-Golay smoothing filter. However, there is no substantial difference in the reference velocity between excitations employing the nonlinear leastsquare alone and those using the Savitzky-Golay smoothing filter. Figure 6 shows that the reference acceleration for all joints is smooth, reducing the shock and effect on the robot arm and ensuring the smooth operation of the ABB irb1600 robot. The black line represents the acceleration after nonlinear least-squares identification with a disturbance, whereas the green line represents the acceleration after the Savitzky-Golay smoothing filter. However, the reference   Ogbemhe et al., Cogent Engineering (2022) acceleration differs significantly between excitations using the nonlinear least-squares method alone and those utilising the Savitzky-Golay smoothing filter. Figure 7 demonstrates that the reference torque for all joints is smooth, eliminating stress and effect on the robot arm and ensuring the ABB irb1600 robot operates smoothly. The black line depicts the acceleration following nonlinear least-squares identification with a disturbance, whereas the brown line represents the torque following the Savitzky-Golay smoothing filter. However, compared to other reference trajectories, the disruption during excitations has the most significant influence on the reference torque. Table 2 shows the Symbolic expressions for the inertia parameters and the estimated values for the dynamic parameters using the Savitzky-Golay smoothing filter.

Conclusions
This study derived the robot dynamic model in its canonical version by employing the Euler-Lagrangian equation. A Freudenstein 1-3-5 trajectory, constructed based on Fourier series expansion, was employed to excite each link of the industrial robot. A slower acceleration rate is attained in the Freudenstein 1-3-5 profiles owing to the higher employment of Fourier series expansion terms. At the same time, limits were enforced to prevent discontinuities in displacement, velocity and acceleration at the start and ends. Nonlinear least-squares based on the Levenberg-Marquart approach was used to find the dynamic parameters of a six-degree-of-freedom ABB irb1600 robot. To boost the identification scheme's performance, we utilise Savitzky-Golay smoothing filters to eliminate a noisy signal with a significant frequency spread. The standard deviation of the torque is compared in Table 3 using nonlinear least squares and Savitzky-Golay filters. The decreased standard deviation values produced with the Savitzky-Golay filters suggest that the values are closer to the mean, which is a desirable attribute.    Ogbemhe et al., Cogent Engineering (2022) m is the mass, and a are link length and centre of mass; f v and f s represents viscous and Coulomb friction parameters respectively and g is the acceleration due to gravity.