Inertial Parameter Identification for Closed-Loop Mechanisms: Adaptation of Linear Regression for Coordinate Partitioning

Modern design and control methods for mechanical systems require accurate simulation models of the system dynamics. Simulation models can be used as a basis for model-based control and estimation, hardware-in-the-loop testing, and verifying control software and design modifications.

In many applications, a rigid multibody model combines sufficient modeling accuracy with high computational efficiency. Rigid multibody models can be built based on kinematic constraints, inertial parameters, and external system forces. The mechanical tolerances of industrially manufactured machines are generally good enough to use kinematic constraints for accurate joint modeling. Inertial parameters, however, are commonly not accurately known, because machines can be equipped with multiple externally purchased components with masses that are not accurately known, or various small components (hoses, guides, plastic covers, etc.) can be added, which are difficult to account for a priori.

Identifying inertial parameters has long been a part of robotics applications, where the established methodology is based on linear regression models [2]. The methodology is based on the fact that robot actuator torques are linear with respect to a specific set of parameters, referred to as standard parameters. This allows the parameters to be identified using a regression model, where joint torques serve as response variables and standard inertial parameters (or their identifiable subset) make up the parameter vector. The regression matrix can be formed based on an inverse dynamic model. The method has been also applied to multibody dynamics framework in Ref. [3], where the method is specifically formulated for Cartesian reference point coordinates.

Because the linear-regression-identification method has been extensively used in robotics [4,5], its requirements and limitations have been the subject of many studies. One of the main challenges of the method in robotics is to reproduce motion data when only encoder measurements are available. The state-of-the-art approach to solving the problem is to filter encoder readings with a zero-phase low-pass filter and produce joint velocities and accelerations via numerical differentiation [4]. In Ref. [6], Fourier-series-based motion trajectories were proposed to enable the analytic computation of joint velocities and accelerations to achieve better signal-to-noise ratios for the computed joint velocities and accelerations.

Motion reconstruction is also studied in multibody dynamics, where the proposed approaches are commonly based on inertial measurement units (IMUs). In Ref. [7], the motion of an excavator is reconstructed using a kinematic Kalman filter and IMUs. Moreover, a recent study [8] presents an approach for motion reconstruction, where it is carried out using IMUs and assuming that the kinematic constraints will be fulfilled.

Vision-based motion-capture techniques are an alternative to IMU-based motion-capture. In Ref. [9], the three-dimensional (3D) spatial motion of the robot end of a three-degree-of-freedom rectangular robot is tracked using a stereo camera system and a visual marker. In Ref. [10], a vision-based motion-tracking technique is presented to track the motion of both a slider-crank and weaving-loom mechanism using camera and visual markers. The approach provides full-field position-level measurement data (i.e., full mechanism motion can be tracked at once) with minimal measurement noise and measurements that can be performed nonintrusively without adding sensors to the mechanism that could modify the dynamics.

The identification methods developed for robotics usually involve specific modeling methods that are hard to use in more general multibody systems. To overcome this inconvenience, [11] proposed a methodology that allows inertial parameter identification while maintaining the constrained DAE structure of general multibody problems. This methodology, however, requires partitioning of the equations of motion into parameter-dependent and parameter-independent parts, which was done in Ref. [11] using symbolic manipulation of the equations of motion.

However, a multibody system, originally described by a system of DAEs, can be transformed into the ordinary differential equation form by applying coordinate-partitioning methods, also referred to as the method based on a projection matrix R [12]. It should be noted that the method in Ref. [12] is agnostic to the coordinate system used and can be applied to various sets of generalized coordinates (e.g., natural coordinates, reference point coordinates, topological joint coordinates, or a mixed set of all of the aforementioned). In this method, coordinates are partitioned into subsets of independent and dependent coordinates. This enables the elimination of the constraint equations from the equations of motion. In simple systems, such as the case example of this study, the independent coordinates can be determined intuitively. For more complex systems, the partitioning can be carried out using Gaussian triangularization or QR or SV decomposition [12].

The objective of this work is to build a conceptual understanding of how the linear-regression-based identification approach can be adopted for closed-loop multibody systems. For generality, the approach is presented using generalized coordinates q and coordinate partitioning to account for motion constraints. This study significantly extends the approach for multibody system identification presented in Ref. [3] by allowing the use of any generalized coordinates, not just Cartesian coordinates. This generalization enhances the applicability and flexibility of the identification process for closed-loop multibody systems. Accordingly, the actual numerical implementation is performed using a double-step semirecursive formulation [13] where joint coordinates are considered as generalized coordinates. According to this scheme, the closed-loop system is temporarily opened by introducing cut-joints [14]. Afterwards, the opened joints are closed again by introducing loop-closure constraint equations in terms of the joint coordinates. Different motion-reconstruction techniques are not covered in detail by this study. However, the method presented in Ref. [10] is used for the experimental validation to produce the kinematic measurements using camera-capture methodology.

The presented approach is validated through both a numerical example and an experimental test of a slider-crank mechanism. The results from the numerical example show good agreement with the original system parameters, and both numerical and experimental results show good agreement between the predicted and measured torque responses of the system. The mechanism used as a validation throughout the paper is depicted in Fig. 1.

where the inertia tensor and first moments of mass are defined with respect to a fixed and known body coordinate system origin. The inertia tensor can be transformed afterwards to the center of mass (X_j, Y_j, Z_j) of a body using the parallel axis theorem.

where H is assumed to be a representative sample of the total energy gradient throughout the state space, r is a number of random coordinates and velocities, and h is the total energy gradient defined in Eq. (2). It is important to note that although the total energy is linear with respect to parameters, it is highly nonlinear with respect to kinematics. Therefore, h, evaluated at different random points in the state space, reveals the overall linear dependencies between the parameters. Since each of the columns of H represent the partial derivative of total energy H with respect to a standard parameter at different points in the state space, the linearly independent columns (nonunique choice) of H determine the inertial parameters required to completely describe the system dynamics.

Some of the columns of H often contain only zeros, which implies that the dynamics of the system are completely independent of the respective standard parameters. For example, the inertia tensor of Body 3 in the Slider-crank mechanism in Fig. 1 does not affect system dynamics, because the body does not have any rotational degrees-of-freedom. These zero columns of H can be eliminated to simplify the identification problem [16]. In this study, the matrix H after eliminating the zero columns is marked as $Hred$⁠.

where $Hred,ind$ and $Hred,dep$ represent independent and dependent columns of $Hred, Kind$ is a vector of independent parameters, $Kdep$ is a vector of dependent parameters, $β$ is a mapping matrix between independent and dependent columns of $Hred$⁠, and $KB$ is a vector of base parameters.

where the dependent parameter vector is set to zeros to form an identification problem with a unique solution (the number of unknown parameters equals the rank of H). However, any other arbitrary choice can be made here, and if other prior information is available, this can be exploited to introduce more representative values. The main advantage of setting the dependent parameters to zeros is that the later introduced partial derivatives with respect to base parameters can be computed as partial derivatives with respect to independent parameters. However, this computational simplification is mathematically equivalent to using the more complex parameter vector, where the term $βKdep$ is also present.

where $Kindold$ is the parameter vector originating from the initial identification, $Kdepnew$ is a vector with desired dependent parameter values, and $Kindnew$ is a vector of updated independent parameter values. This update procedure can be also used to transform the solution to a minimum norm form, where $||K||$ is minimized.

where R is a velocity transformation matrix, Q is a vector of generalized forces without the actuation forces, $Qact$ is a vector of the actuation forces, M is a mass matrix, and $z˙˙$ is a vector of independent accelerations. The distinction between $Qact$ and Q is related to their parameter dependency and measurability. The externally applied forces that are measurable, are collected to the vector $Qact$ while the other generalized forces, dependent on parameters to be identified, are present in Q. For scleronomic constraints, $Sc=R˙z˙$ represents the generalized accelerations with zero independent accelerations [12], and $Φ$ is a regression matrix.

In Eq. (11), the system must be actuated to achieve persistent excitation of the dynamics to enable accurate identification. In other words, the identification scheme cannot identify the parameters of a system affected by gravity alone, because an identification problem where $Qact=0$ has the trivial solution of $KB=0$⁠. However, if some of the base parameters are a priori known, they can be moved to the other side of Eq. (11) together with the respective columns of the regression matrix $Φ$⁠. This way, an identifiable problem may also be formulated in the case of a purely gravity-actuated system.

The partial derivatives of the mass matrix and generalized force vector with respect to standard parameters depend on the types of coordinates used, and the expressions should therefore be derived separately for different types of generalized coordinates. However, the derivation simplifies, because the end result is known to be independent of the inertial parameters, which allows all remaining inertial parameters to be set to zero after differentiation. This procedure may be necessary when computing partial derivatives of different matrix elements with respect to the mass of a body m_j, which can result in remaining matrix elements with centers of mass X_j, Y_j, and Z_j.

The identification algorithm is depicted in Fig. 2. As the figure shows, the algorithm assumes that the independent coordinates and actuation torques are measurable and uses them as inputs. Because numerical differentiation has a tendency to amplify measurement noise levels, the measured independent coordinate data are filtered before numerical differentiation is carried out. Because the identified system is modeled using a rigid-body model with no inequality constraints, higher frequencies, typically originating from deformation and impacts, can be considered as noise. Filtering in this study is applied using a zero-phase low-pass filter with Butterworth tuning strategy and is implemented using the Matlab filtfilt function [19], as in Ref. [4]. Also as depicted in the figure, filtering and numerical differentiation can be considered preprocessing phases of the algorithm, because they are purely signal processing techniques that do not require any information about the system.

The inverse dynamics model (Eq. (11)) is formed after the independent coordinates, velocities, and accelerations of the system have been evaluated. In this process, the set of independent coordinates z are first transformed to a complete set of coordinates q by solving the position problem of $C(q)=0$⁠. C is the set of constraint equations. Second, the independent velocities are transformed to a set of generalized velocities using the velocity transformation matrix R. The R-matrix, in turn, is computed from the Jacobian of constraints (C) with respect to generalized coordinates. Its computation thus requires the previously computed vector of generalized coordinates. Once these quantities are obtained, the remaining matrices depending on the system's kinematic state can be computed $(R˙z˙$⁠, $∂M∂KB$⁠, and $∂Q∂KB)$⁠.

These matrices, together with the acceleration and actuation force data, form the linear regression model for a single time-step. In the last step of the algorithm, the individual regression models from multiple time instances are combined, and the parameter set for that combined model is found using an ordinary least squares approach.

The proposed identification algorithm was validated using a numerical example for the slider-crank mechanism shown in Fig. 1. The body-fixed coordinate systems were attached to the revolute joint locations. The mechanism was modeled using the relative joint coordinates: q₁ between the ground and Body 1, q₂ between Body 1 and Body 2, and q₃ between the ground and Body 3. A cut-joint was applied for the revolute joint between bodies 2 and 3. Because the mechanism has one degree-of-freedom, one independent coordinate must be defined. In this study, the chosen independent coordinate is z₁ = q₁.

The parameters of the system are shown in Table 1. In the table, the independent parameters are highlighted dark gray and the dependent parameters are highlighted light gray. The nonhighlighted parameters do not affect system dynamics.

For the numerical validation, the data were produced by simulating the mechanism using Simscape Multibody [20]. The mechanism was kinematically driven, with a kinematic reference set for crank angle as crank torque was computed. The crank angle together with crank torque was used as artificial measurements augmented with artificial measurement noise $σz=1×10−3$ rad, $σT=1×10−2N·m$⁠, respectively, and input to the identification algorithm. The reference crank angle (q₁ precise), noisy crank angle (q₁ “measured”), and filtered crank angle (q₁ filtered) during the simulation are depicted in Fig. 3. Filtering of the crank angle data was performed using a zero-phase filter with a Butterworth tuning and a passband frequency of 10 Hz and a stopband frequency of 20 Hz. These frequencies were chosen to be substantially higher than the excitation frequency. The left half of the figure (showing training trajectory) depicts the encoder data used to find the parameter set that best fulfills the regression model (Eq. (11)). Correspondingly, the right half of the figure (showing validation trajectory) was used to test the model's ability to predict crank torque based on measured crank angle. Figure 4 shows the results of the measured and predicted actuation torques together with the original noiseless torque. Similarly as in Fig. 3, left half data were used to train the regression model, while the rest of the data were used for validation. As the figure shows, even with noisy torque measurements, the regression model can identify a model that is capable of predicting the original joint torques with high accuracy.

Table 2 presents the reference and identified values of the base parameters. The reference parameter values are computed by substituting the parameter values used in simulation into Eq. (9). As the table reveals, the identified parameter values show good agreement with the real reference values. However, the parameter values in the table do not correspond to the values given in Table 1 because of the lumping of dependent parameter values into independent parameters, as indicated in Eq. (9).

To further validate the practical feasibility of the introduced algorithm, the proposed method was tested using the experimental setup shown in Fig. 5. The motion of the system was recorded using the camera-capture technique reported in Ref. [10]. The movement of the mechanism was captured at a frame rate of 800 fps using a stereo camera setup comprising a JAI SP-12000M-CXP4 camera with ZEISS Distagon T* 2/28 mm (industrial edition) lens in combination with a Ximea xiB-64 CB120RG-CM-X8G3 camera with ZEISS Distagon T* 2/28 mm lens. The stereo camera setup ensured the extraction of 3D spatial motion by camera triangulation [21].

To facilitate the motion tracking, visual fiducial markers were applied at the revolute joint locations of the mechanism together with an ArUco marker on Body 2 as shown in Fig. 5. These fiducial markers are subsequently tracked in all the frames of the camera recording by exploiting the coarse-to-fine tracking approach as reported in Ref. [10].

First, the motion of the mechanism at each point in time was estimated using an affine geometric projection of the initial system configuration. The same matrix needed to project the detected ArUco marker locations [22] from the initial (reference) time instance to the current time instance was imposed. Then, this coarse motion estimate was refined with subpixel accurate image correlation techniques [23] and triangulated to the result in the 3D position-level measurement data. By using the initial system configuration as a fixed reference, drifting errors in the motion tracking were avoided.

where F_c is the Coulomb friction coefficient. In this work, the friction dependency on normal forces, or more generally on reaction forces, is deliberately neglected because that dependency would introduce parametric nonlinearities for the regression model preventing the use of linear regression for parameter identification.

In the experimental verification process, the mechanism's model was identified using data from one experiment, while data from a separate experiment served to verify the model. Practically, identification involved finding a model that most accurately predicts the measured torque response, based on the observed kinematics. Conversely, verification entailed comparing the measured torque from the dataset of the separate experiment against the torque predictions made by the identified model, again based on kinematic measurements. Notably, in the experimental identification phase, prior knowledge about the parameter values was applied. This was achieved by setting the values of $m1Y1$ and $m2Y2$ to zero, acknowledging the symmetry of Body1 and Body2 with respect to the X-axis.

Figure 6 illustrates the motion trajectory used in the experimental setup for training the regression model. In the figure, the mechanism initially accelerates from rest at a constant rate, followed by a deceleration phase that mirrors the initial acceleration. This process is designed to activate the acceleration dependent inertial parameters (such as $I1,zz$⁠) of the system while the variation in velocities aims to develop a model that is effective across different velocities. Although the trajectory shown in the figure appears smooth due to the original camera measurements, further filtering is necessary to avoid high noise levels in the computed acceleration data. The filter was tuned using Butterworth tuning strategy with bandpass frequency of 50 Hz and stopband frequency of 60 Hz. Filter was tuned by increasing the bandpass frequency until degradation was observed in identification results, which was monitored detecting the value of error variance (error being difference between measured and predicted torque).

Figure 7 depicts the system's torque response from the experiment used to identify the model together with the response predicted by the identified model. As the figure illustrates, a clear discrepancy between predicted and measured torques is present at the startup and stop, when the mechanism operates at low angular velocity. A probable reason for this is the complex nature of friction forces at low velocities, which cannot be modeled accurately using the simple Coulomb friction model. With higher rotational speeds, the correspondence between measured and predicted torques can be considered acceptable. To verify the identified model, it was applied to a separate experimentally collected dataset. Figure 8 shows motion data of the mechanism used to verify the identification results. In this verification experiment, the mechanism was accelerated rapidly to speed, which was then held approximately constant for the remainder of the experiment. Figure 9 depicts the torque measured from the verification experiment together with the torque predicted using the motion data from the same experiment and the parameters identified from the training experiment. The figure demonstrates good agreement between the measured and predicted torque, with the exception of the startup phase.

The parameters identified from the training dataset (Fig. 7) are presented in the Table 3. For the experimental setup, the exact parameter values were a priori unknown and were supposed to be different from the parameters used in the numerical example. However, the identified parameter values presented in the table seem reasonable for the given experimental setup. It should be noted that the parameters $m1X1$ and $m2X2$ were not identified and were set to zeros utilizing a priori information of mechanism's geometrical symmetries.

The aim of this work is to bring regression-model-based parameter identification methodology into a generalized multibody framework, where the motion of bodies is restricted by various kinds of constraints. With over three decades of research, the methodology is well established in the field of robotics.

The results show that the method is capable of accurately identifying the underlying rigid-body model parameters of a slider-crank multibody system when using simulation-based measurement data. Moreover, the experimental results from a similar setup show that the identified model is capable of predicting the actuation torque with acceptable accuracy. In addition to the inertial parameters, the experimental setup also included the identification of the Coulomb friction coefficient of the slider.

In addition to bodies, linear springs and dampers are elements commonly incorporated in multibody systems. Integration of these elements to the presented identification framework is straightforward, as demonstrated in Ref. [3], because like standard inertial parameters, the relationship between force and spring/damping coefficient is linear.

The presented method eliminates and groups some of the inertial parameters to achieve an identifiable problem. It is important to note that this simplification is a requirement for the given input–output data because the parameters having linearly dependent effects on the system's response are fundamentally unidentifiable. However, this limitation makes the approach incapable of providing reaction forces or properties of the individual bodies. These limitations make the approach more suitable for state estimation or control applications rather than detailed design analysis.

In future work, the camera capture methodology could be paired with measurements from inertial measurement units (IMUs) using information fusion techniques such as Kalman filtering. By changing the filtering strategy from the current noncausal zero-phase low-pass filter to a causal Kalman filter, parameter identification could be possibly executed online, shifting away from the current offline scheme. Moreover, the method could be utilized for larger scale systems, such as mobile machinery, where the effect of inertial parameters to system's force response is more significant.

The first author would like to acknowledge KU Leuven Mechatronic System Dynamics (LMSD) and in particular Professor Frank Naets for hosting his research visit during which the main ideas of this study were developed.

• Flanders (FWO) (Grant No. G095120N; Funder ID: 10.13039/501100003130).

• Thijs Willems (Grant No. V462123N; Funder ID: 10.13039/501100003130).

• Flanders Make, the strategic research centre (Funder ID: 10.13039/100020468).

• KU Leuven (Funder ID: 10.13039/501100004040).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.