The dynamic model of a serial manipulator, relating robot motion to joint driving torques, plays a crucial role in designing advanced control laws. For an n-DoF serial manipulator composed of rigid links, the dynamic equation can be derived through first principles, such as the Lagrangian method. Ignoring external disturbances and joint frictions, the dynamic model can be generally expressed in the following form:where are the vectors of generalized coordinates, generalized velocities, generalized accelerations, respectively; M, C, respectively represent the mass, the Coriolis and centrifugal matrices; G represents the gravity term; is the vector of generalized torques/forces. For the n-DoF serial manipulator, the process of first-principles modeling is tedious and skillful, and inevitably involves ideal assumptions. In contrast, SINDy-based model reconstruction formulizes the dynamic behavior of the manipulator using data, avoiding the ideal assumptions as well as mathematical derivation and parameter identification in modeling process.

The model of the n-DoF serial manipulator has the following property.

Property 1. The left side of Eq. (1) is linear with respect to system parameters [43], i.e.,where ψ represents the regression matrix and θ the vector of system parameters. The fundamental principle of SINDy is to sparsely represent a function as a linear combination of nonlinear terms. To this end, ψ is to be replaced by the libraries of candidate functions, namelywhere is the library of candidate functions, is the sparse matrix of coefficients. According to the characteristics of rotational joints and joint frictions of the serial manipulator, consists of constants, trigonometric functions, exponential functions, and polynomials. Subsequently, by recording the signals of robot motion and joint torques, we obtainwhere is the sparse vector of coefficients of ith joint, p the number of candidate functions, and i = 1, 2, ..., n. At last, is determined by solving Eq. (4) in different regression strategies.

As the DoF of the serial manipulator increases, more candidate functions with complicated form will appear in the library and the size of will grow exponentially, eventually leading to curse of dimensionality especially for the case where q, , and all have significant measured values. Furthermore, signals of higher time-derivative measurement are more vulnerable regarding the signal-to-noise ratios (SNRs), implying that the regression may be severely disturbed by the noise if the signals of q, , and are simultaneously taken into account for regression. In other words, separately collecting data according to different working conditions may facilitate efficient model reconstruction.

Considering typical tasks of serial manipulators, stages with different data features will appear alternatively, indicating the possibility of quantifying such data features according to the activation of different signals. To this end, we first define the RAI as follows:where a ∈ can be substituted by the vector of generalized velocities or generalized accelerations , ‖·‖₁ is the sum of the magnitudes of the vector in space, and W is the weighting coefficient of the activation level. S(·) is the function of activation switch that takes the formwhere h is the coefficient that adjusts the sensitivity of the activation switch. After processing the generalized velocity and generalized acceleration signals with RAI, the collected data can be clustered using the k-means method. The procedure of the data clustering based on RAI is presented in Algorithm 1.

Figure 1 illustrates the process of data clustering based on RAI as sketched in Algorithm 1. Firstly, the n-DoF serial manipulator performs a sorting operation that involves typical states of the joints such as static, uniform velocity, and variable velocity. To simulate the noise in measurement, we add Gaussian White noise with SNRs of 60 and 50 dB to the signals of generalized velocities and generalized accelerations, respectively. Signals from each joint are collected during the operation of the robotic manipulator. Subsequently, the collected data is pre-processed by calculating RAI for the vector of generalized velocities and generalized accelerations andwhere the parameters are chosen as h = 60, W₁ = 0.6, and W₂ = 0.32. Finally, the k-means method is employed to cluster the processed data as detailed in Algorithm 1. It’s worth noting that the number of clusters k is selected by the elbow method. Here the within-cluster sum-of-squares (WCSS) is introduced as an evaluation index of clustering, which is defined as [44]where ‖·‖₂ is the distance of the vector coordinate from the origin of the vector space. As shown in Fig. 1, the elbow point of the number of clusters arises when k = 3. Additionally, Algorithm 1 is run 20 times with randomly initialized clustering centers, and the maximum number of iterations is 12 while the average number of iterations is 5.75, indicating the effectiveness of the data clustering.

Remark. The reason that RAI defined by Eq. (5) is termed as “relative” is as follows. The first term on the right-hand side of Eq. (5), i.e., S (a), appears as a filter that determines which signal is considered as activated, and the second term quantifies the activation level based on which the post-clustering data is uniformly distributed around S (a). Therefore, Eq. (5) indicates the relative magnitude of the activation. However, it should be noted that the first term in RAI is still dominant if a proper value is assigned to W.

In this subsection, we aim to elucidate the physical interpretation of the clustering results observed in the last subsection. Backtracking from three clusters to collected data, they respectively represent the static state, constant velocity motion, and variable velocity motion as shown in Fig. 2. Data in the first cluster leads to the least significant RAI and depends almost only on generalized coordinates, implying that the library of candidate functions that reconstructs the model for the static case has the simplest form and the smallest size. Meanwhile, data in the second and third clusters activate respectively the functions involving generalized velocities and generalized accelerations, indicating that the complexity of candidate functions will be growing in the corresponding cases.

Revisiting the dynamic equation of the serial manipulator and ignoring the joint frictions, the motion signals in the three cases, namely, the static state, motion with constant velocity, and motion with variable velocity, will sequentially activate the gravity terms, the Coriolis and centrifugal terms, and the inertia terms, respectively. Consequently, reconstructing the dynamic model of the serial manipulator in a stepwise manner, based on the recognition of the aforementioned cases, can effectively reduce the complexity of symbolic regression for (i) the regression results of the current step can be inherited by the subsequent step, and (ii) the complexity of matrix operations, essentially abundant in the SINDy regression, can be significantly reduced. Furthermore, since the signals of generalized velocity and generalized acceleration are separately used in different steps, this approach can effectively mitigate the impact of noise on the reconstruction results.

Combined with the domain knowledge, the library of candidate functions i.e., is assumed to be constructed by the elementary functions including generalized velocities, generalized accelerations, generalized coordinates, and trigonometric functions of generalized coordinates as building blocks. Due to the kinematic constraints inherent in multi-DoF serial manipulators, expressions concerning generalized coordinates in the dynamic model typically exhibit greater complexity when contrasted with expressions in terms of generalized velocity and generalized acceleration. To fully describe the constituents of library and reduce non-essential candidate functions, DLG is proposed. For the library generation of SINDy, the first feature set is given as

In Eq. (10), c_1i is the elementary function, the basic components for feature set , andwhere E_i = {sin q_i, cos q_i} if the ith joint is a revolute joint and E_i = {q_i} if the ith joint is a prismatic joint. N₁ represents truncated order of the elementary functions, which determine the complexity of features in Since the kinematic relations are simpler in the latter, prismatic joints are more friendly to library generation. {1} in Eq. (11) is to ensure that as N₁ increases, the feature set still contains some simple terms. For instance, if N₁ = 3 in Eq. (10), the complex features in may be {sin q₁cos²q₂} and the simplest feature is {1}. We definewhere |·|, called cardinality of the set, denotes the number of elements of the set. The number of features (i.e., elements)in can be calculated as

Then, the second feature is given aswhere

We define

The number of features in can be calculated as

Finally, the complete feature set then has the following form:

Each term on the right-hand side of Eq. (2) can be represented by a linear combination of the features (i.e., elements)in , therefore all features in make up the library of candidate functions. The procedure of DLG is presented in Algorithm 2. DLG excludes non-essential candidate functions containing higher-order generalized velocities and higher-order generalized accelerations, which is beneficial for both SINDy and the stepwise SINDy.

Inspired by the clustering results together with the terms appearing on the left-hand side of Eq. (2), we conclude that the model of the n-DoF serial manipulator has the following property.

Property 2. The right-hand side of Eq. (2) can be rewritten as

With this, the dynamic model of the n-DoF serial manipulator can be reconstructed in three steps. Herein, we denote the symbols in the ith step by {·}^S_i with i = 1, 2, 3. In SINDy-based model reconstruction, ψ_i is replaced by the libraries of candidate functions, namely

When performing DLG in stepwise SINDy, elementary functions of the first feature set is the same as Eq. (11) and we have

Meanwhile, elementary functions of the second feature set require a specific design for each step.

(1) Step I

In the first cluster, data with the least significant RAI depends almost only on generalized coordinates, implying that the library of candidate functions that reconstructs the model for the static case has the simplest form and the smallest size. Therefore, we choose the data collected from the static case for the first step of the model reconstruction. Then we have

In this step, we use the following replacement:

Therefore, the elementary function of the second feature set in Step I is selected as and the second feature set is determined as

(2) Step II

In the second cluster, the data activates functions involving generalized velocities. The data of motion with uniform velocities is collected for the second step of the model reconstruction. Then it is clear that

where Ψ₁ (q) θ₁ is inherited from the regression results of the first step. In this step, we use the following replacement:with a particular form of the second feature set given aswhere

(3) Step III

In the third cluster, the data further activates the functions involving generalized accelerations. The data of motion with variable velocities is collected for the final step of the model reconstruction. It is clear thatwhere Ψ₁ (q) θ₁ and are inherited from the regression results of the first and second steps, respectively. In this step, we use the following replacement:with a particular form of the second feature set given aswhere

In this subsection, the complete procedure of model reconstruction for the n-DoF serial manipulator is presented, and the flow diagram is shown in Fig. 3. At first, signals conforming to the three different clustering patterns are collected separately after determining the range and number of data collections, and the strategies for the data collections are proposed. In the data collection of the third clustering pattern, the form of Fourier-series-based trajectories is given aswhere i = 1, 2, ..., n, k and ω are randomly generated. Next, the optimal truncated orders of feature sets need to be determined. The guiding principle for selecting the optimal truncation order is to minimize it while ensuring an adequate description of the dynamic behavior of the system, and we will give examples in the next section. Then, libraries of SINDy and stepwise SINDy are generated by using DLG. Finally, the reconstructed models of SINDy and stepwise SINDy are obtained by solving the regression problem. Herein, the Lasso regression with L1 norm penalties/constraints [45] is employed for finding sparse solutions to regression problems. The problem is formulated aswhere λ is the sparsity promoting parameter. For stepwise model reconstruction, the reconstructed model obtained in the current step needs to be passed along to the next step, and the reconstructed model in the ith step is given by

SCARA is a prevalent industrial robot system that is widely recognized for its high rigidity and positioning accuracy in the horizontal plane [46]. The SCARA robot features a kinematic structure where the first two revolute joints have parallel axes, enabling planar motion of the end effector in the horizontal plane. Additionally, the vertical position and orientation of the end effector are independently controlled by the third prismatic joint and the fourth revolute joint, resulting in minimal kinematic coupling among the four DoFs and a simplified dynamic model. In this subsection, a typical 4-DoF SCARA manipulator with three revolute joints and one prismatic joint is presented as the first example to test the performance of the proposed stepwise model reconstruction approach. The configuration and parameters of the SCARA manipulator are shown in Fig. 4 where the simulations of inverse dynamics are performed in MATLAB-Robotic Tool-box.

The dynamic model of the SCARA manipulator is reconstructed following the procedure shown in Fig. 3. At first, the motion range of data collections is prescribed as shown in Table 1, and the numbers of data collected for the stepwise SINDy across three steps are 2000, 4000, and 4000, respectively. Then, for the SCARA manipulator with a simple dynamic model, the optimal truncated orders of the feature sets are small and can be selected in accordance with the guiding principle. The comparison between SINDy and stepwise SINDy of the SCARA manipulator is shown in Table 2. Under the optimal truncation orders, the total number of candidate functions for stepwise regression is 168, and for SINDy is 360, representing a decrease of 53.33%. The reduction in the size of the library of candidate functions is accompanied by an improvement in computational efficiency. Performing a single calculation of least squares regression, the total time consumption for the stepwise SINDy is 0.1165 s, whereas for SINDy it is 1.48 s, representing a decrease of 92.13%.

In addition, we study the performance of the stepwise SINDy in the presence of data noise. To simulate the noise in measurement, Gaussian White noise with SNRs of 60, 40, and 20 dB are added to generalized coordinates, generalized velocities, and generalized accelerations, respectively. The sparse vectors of coefficients, i.e., i = 1, 2, 3, 4, are solved from Eq. (32) where λ is chosen to be 0.02. Components in with an absolute value exceeding 0.001 are considered as the active items. As shown in Fig. 5, the number of active items in the reconstructed model of the stepwise SINDy is significantly reduced compared to that of SINDy, by 87.1% for , 64.3% for , 88.6% for , and 88.2% for . The fewer number of active items indicates that the reconstructed model is more concise, which generally makes the model more predictive. To compare the predictive ability of the reconstructed models of the two methods, a trajectory outside the range of data collection is used for torque/force predictions. As shown in Fig. 6, the reconstructed model using the stepwise SINDy provides a better predictive ability than the reconstructed model using SINDy.

To conveniently compare with the experimental results obtained from a serial manipulator with three revolute joints, a three-axis serial manipulator is introduced to validate the effectiveness of the stepwise model reconstruction. The diagram and parameters of the three-axis serial manipulator are shown in Fig. 7 where the simulations of inverse dynamics are performed in MATLAB-Robotic Toolbox.

Similarly, the dynamic model of the three-axis manipulator is reconstructed following the procedure shown in Fig. 3. The motion range of data collections is presented in Table 3 and the numbers of data points collected for the stepwise SINDy across three steps are 1000, 10000, and 10000, respectively. Despite having only three degrees of freedom, the three-axis serial manipulator involves multidimensional motions, thus leading to a dynamic model more complicated than that of the SCARA manipulator. To find the optimal truncated orders, we generate libraries of candidate functions for SINDy and the stepwise SINDy at different truncated orders, and then calculations of least squares regression for the collected data sets are performed. The regression errors at different orders are shown in Fig. 8. The figures show that the errors of least squares regression decrease rapidly when the orders reach a certain value, and no longer decrease remarkably when the orders keep increasing, implying the existence of optimal truncated orders. The comparison between SINDy and stepwise SINDy of the three-axis manipulator is shown in Fig. 9. Under the optimal truncation orders, the total number of candidate functions for stepwise regression is 2968 and SINDy is 5880, representing a decrease of 49.52%. Performing a single calculation of least squares regression, the total time consumption for the stepwise SINDy is 32.3 s, whereas for SINDy it is 499.9 s, representing a decrease of 93.54%.

Furthermore, we discuss the performance of the stepwise SINDy under data noise through the following three cases. In Case I, all collected data is noise-free. In Case II, we add Gaussian White noise with SNRs of 70, 50, and 30 dB to the signals of generalized coordinates, generalized velocities, and generalized accelerations, respectively. In Case III, the collected data is added with stronger Gaussian White noise with SNRs of 60, 40, and 20 dB to the signals of generalized coordinates, generalized velocities, and generalized accelerations, respectively. The sparse vectors of coefficients, i.e., , i = 1, 2, 3, are solved from Eq. (32) where λ is chosen to be 1. Components in with an absolute value exceeding 0.01 are considered as the active items. As shown in Fig. 10, the number of active items in the regression results of the stepwise SINDy is close to that of the traditional SINDy in Case I, indicating an equal level of sparsity of applying the two methods in the absence of noise. In Case II, the number of active items in the regression results of the stepwise SINDy is significantly reduced compared to that of the traditional SINDy, by 75.62% for , 83.81% for , and 50.00% for . With the enhancement of the noise as in Case III, the sparsity of regression results of stepwise SINDy turns out to be more pronounced as reductions of 83.34% for , 85.57% for , and 85.78% for are observed. A reduced number of active items in the reconstructed model signifies a higher level of sparsity, thereby typically enhancing its predictive capability. To compare the predictive ability of the reconstructed models of the two methods, a trajectory outside the range of data collection is used for torque predictions. As shown in Fig. 11, the reconstructed model using the stepwise SINDy still provides promising prediction in Cases II and III, while the reconstructed model using the traditional SINDy fails to accurately predict the torques in Case III with the strong noise.

As shown in Fig. 12, an experimental platform of three-axis serial manipulator is built to further validate the effectiveness of the proposed method under realistic working conditions. The manipulator is equipped with three Elephant^@ joint modules, each of which is integrally composed of a servopack, encoder, and reducer. The produce numbers of joint modules are MS32 for joint 1 and MS20 for joints 2 and 3. The SpeedGoat^@ Unit is implemented as the programmable controller. In the experiment, we collect data for model reconstruction and verification in the speed mode of the joint modules. Based on the reconstructed model, we design the control law which is experimentally deployed in the torque mode to guarantee the performance of trajectory tracking. The sampling frequency of joint angles and joint velocities is 1000 Hz.

In the experiment, one has to consider the friction in the joint module in which the reducer may cause complicated frictional behaviors. In addition to the aforementioned procedure of data-driven model reconstruction (as shown in Fig. 3), the friction is modeled by the typical Stribeck model [47], given by

Note that the sign function and exponential function in Eq. (34) are difficult to handle in constructing the library of candidate functions. For this reason, we include in the model reconstruction the following treatments: (i) for the sign function, we categorize the collected data into the forward stage and backward stage according to the direction of the joint motions and then perform the stepwise SINDy for either stage of the data, and (ii) for the exponential function, we add exponential functions of joint velocities with different coefficients, i.e., where i = 1, 2, 3 and v_s = 0.01, 0.02, ..., 1, to the library of candidate functions. It is worth mentioning that even in the static state the friction may also incorporate with the gravity, resulting in the difficulty of experimentally acquiring effective data for Step I in the stepwise SINDy. Alternatively, in this step the joint modules are actuated at a very low constant speed, i.e., 0.1 °/s, to overcome the coulomb dry friction. As a complement to Fig. 3, Fig. 13 illustrates the procedures of data collection and model reconstruction for the first two steps of stepwise SINDy in the experiment.

The reconstructed models are employed to predict torques and the results are shown in Fig. 14. The reconstructed model in the first step captures the characteristics of the gravity and the joint dry frictions of the mechanism, so the torque predictions for the first and third joints, being less affected by gravity, are closer to the dry friction. The model reconstructed in the second step further predicts the velocity-dependent dynamics of the mechanism, such as the velocity terms in the friction and the Coriolis and centrifugal terms. The last step reconstructs the complete dynamic model of the three-axis serial manipulator. The experimental results demonstrate that as the number of steps increases, the reconstructed model obtained at each step of the stepwise SINDy becomes more accurate in predicting the joint torques. Additionally, we design the corresponding control strategies based on the reconstructed models obtained from Step I and Step III of stepwise SINDy, respectively given bywhere e is the vector of tracking error of joint angles, and is calculated from the reference signals (RS). The controller with the model reconstructed in Step III is designed aswhere , and are also calculated based on the RS. To demonstrate the effectiveness of the proposed model reconstruction method, controllers Eqs. (35) and (36) are to be compared with the classic PID control in the torque mode. The gains of all controllers are identically chosen as K_p = 80, K_i = 8, and K_d = 2. The external trajectory tracking performance of the three controllers is shown in Fig. 15 where the time axis does not represent continuous time. The interval 1-20 s represents stable trajectory tracking under PID control, 20-40 s under the controller in Eq. (35), and 40-60 s under the controller in Eq. (36). The results of the first period of trajectory tracking show that the PID controller has difficulty in overcoming the dry friction during the joint velocity commutation, leading to the overall tracking performance degradation. Denote by the average of absolute tracking errors of three joint angles with unit of deg. Compared to the model-free PID controller (η₁ = 1.10, η₂ = 1.21, η₃ = 1.23), the controller that uses the reconstructed model in Step I of stepwise SINDy to generate the dynamic compensation significantly improves the trajectory tracking accuracy (η₁ = 0.73, η₂ = 0.50, η₃ = 0.81), which is further improved after implementing the reconstructed model obtained from Step III of stepwise SINDy in forming the dynamic compensation (η₁ = 0.25, η₂ = 0.25, and η₃ = 0.27). It is worth noting that the interpretability of the reconstructed models in the stepwise SINDy supports the design of advanced model-based control laws.