{'fraction_non_alphanumeric': 0.0336080659358246, 'fraction_numerical': 0.00712170921021045, 'mean_word_length': 4.50418410041841, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simplification of multibody models by parameter reduction

Model selection methods are used in different scientific contexts to represent a characteristic data set in terms of a reduced number of parameters. Apparently, these methods have not found their way into the literature on multibody systems dynamics. Multibody models can be considered parametric models in terms of their dynamic parameters, and model selection techniques can then be used to express these models in terms of a reduced number of parameters. These parameter-reduced models are expected to have a smaller computational complexity than the original one and still preserve the desired level of accuracy. They are also known to be good candidates for parameter estimation purposes. In this work, simulations of the actual model are used to define a data set that is representative of the system's standard working conditions. A parameter-reduced model is chosen and its parameter values estimated so that they minimize the prediction error on these data. To that end, model selection heuristics and normalized error measures are proposed. Using this methodology, two multibody systems with very different characteristic mobility are analyzed. Highly considerable reductions in the number of parameters and computational cost are obtained without compromising the accuracy of the reduced model too much. As an additional result, a generalization of the base parameter concept to the context of parameter-reduced models is proposed.


Introduction
The requirement of fast simulations to speed up design or to enable real-time applications has been one of the main concerns of researchers in Multibody Systems Dynamics (MSD) over recent decades. This, in turn, has made model reduction a fundamental topic in the literature on MSD. Here, model reduction most frequently means model order reduction (MOR), that is, a reduction in the number of degrees of freedom (d.o.f.) of the model without compromising its accuracy too much. A common situation is that of linear model order reduction of the elastic d.o.f. of bodies. The most popular techniques include: Component Mode Synthesis, Moment Matching Krylov subspaces, SVDbased and Proper Orthogonal Decomposition-based techniques, etc. For an overview of these techniques see [1]. MOR techniques can also be applied to full multibody systems or subsystems using, for example, the Global Modal Parametrization [2].
Real-time performance is a requirement in many multibody applications, such as feed forward and predictive control [3,4], hybrid (analytical and data-based) model-based control [5], model-based nonlinear filtering [6], modelbased condition monitoring [7,8], etc. Frequently, models do not require body-level flexibility to be accounted for and, therefore, MOR is not suitable. However, due to the strict demands of real-time performance, researchers keep striving to minimize the computational demand required by the models. In this context, formulations that solve the exact dynamic equations using a minimum number of operations as recursive O(n 3 ) and O(n) formulations play a predominant role [9].
Symbolic multibody methods can be used to implement recursive formalisms. This reduces the operation count even more due to the exploitation of symmetries [10]. In addition to this, the development of special inertial parametrizations such as the barycentric inertial parameters [11,10] and minimum or base inertial parameters [12] should be mentioned, as these are used to further lower the operation count. Some heuristics have also been proposed in order to determine more efficient kinematic parametrizations [13].
Other methods further lower the operation count at the expense of accuracy. They remove or approximate some terms of the model functions, or change their computation frequency. These procedures generally rely on symbolic methods. In [14] the authors propose to remove the terms related to Coriolis and centrifugal forces and to take out the kinematics of the terminal part from the dynamic computations; in [15] different terms of the dynamic model are updated using different frequencies; in [16] a statistical analysis of the values of different elements of the dynamic model functions is carried out in order to replace some of them with suitable constants; and in [17] the authors use a spline approximation to determine the dependent coordinates as a function of the independent ones for some closed-loop subsystems (car suspensions). Obviously, some of these strategies can be used simultaneously.
Due to the linearity of the equations of motion with respect to the inertial parameters [18], the typical inertia parameter estimation problem [19] takes the form of a standard linear regression problem where ϕ is the vector of inertial variables, χ is the vector of external forces at different measurement points, and W is the observation matrix that is a function of the extended state at the different measurement points. Standard linear model selection techniques can be regarded as methods that "find the subset of parameters (ϕ) that best fit a given data vector (χ) as a linear combination of some regressors (columns of W)". This selection is a combination of two elements: 1) an error criterion to measure the accuracy of the fit [20], and 2) a heuristic to select the most promising parameter combinations [21,22,23,24,25]. Validation techniques which measure the error criterion for the selected model on a validation data set are frequently used.
In this paper, we propose the use of model selection techniques to obtain approximate parameter-reduced dynamic models of multibody systems. The regression problem is fed with data obtained from simulations based on the exact model and model selection techniques are used to choose the simplest set of parameters that fit these data for a given error tolerance.
To that end, we propose non-dimensional normalized error measures to determine the performance of direct and inverse parameter-reduced dynamic models. Also, as a blind search would result in an enormous computational effort, three different heuristics -QR, Backward Elimination and Forward Selection-are proposed for the selection of the candidate parameter sets. The performance of the proposed methods is tested on two examples of rigid multibody systems. As the high/low mobility of bodies -defined as their ability to exhibit large/small excursions in their respective motion space-plays an important role in the relevance of different dynamic parameters, two multibody systems with very different characteristic mobility are analyzed. Trajectory optimization is used to define sufficiently exciting characteristic data sets for the regressions. Model errors and their algebraic operation count are reported for reduced models with a varying number of parameters obtained with the different heuristics proposed. Based on these results, we discuss the relevance of the proposed methods in the contexts of dynamic parameter estimation and computational complexity reduction in MSD.
The article is organized as follows: Section 2 presents the methods and algorithms proposed in this article. First, in Section 2.1, the classical base parameter determination [12] algorithm is proposed as a parameter reduction method. A brief explanation of the algorithm is given, emphasizing this perspective. In Section 2.2, the parameter reduction method proposed in this paper is introduced as a natural extension of this view. The proposed heuristics are described in sections 2.3, 2.4 and 2.5. In Section 3, the algorithms introduced in this paper are applied to the selected examples. A detailed description of the analyzed systems and the results obtained are presented in Sections 3.1 and 3.2. Section 4 analyzes the results obtained, putting the relevance of the methods proposed into perspective. Possible improvements and a more adequate normalized error criterion for direct dynamic models are proposed and discussed. Finally, in section 5, the major outcomes of this work are highlighted.

Parameter Reduction Method
Using a set of independent coordinates z, the Lagrange equations of a multibody system can be arranged as where d z is the so-called inverse dynamics model (IDM). The IDM makes it possible to obtain the vector of externally applied generalized forces τ z , given the extended state triad (z,ż,z) and the vector of dynamic parameters ϕ of the system. In the context of rigid-body dynamics, symbolic implementations of recursive O(n 3 ) formulations constitute the state-of-the-art general purpose techniques for the efficient computation of the IDM. The direct dynamics model (DDM) is usually obtained by splitting out d z as where M zz represents the mass matrix function and δ z the generalized Coriolis, centrifugal and constitutive forces function. Note that based on these DDM functions, given the state (z,ż) of the system and the external generalized force vector τ z , the generalized accelerations can be obtained as follows: Regarding computational efficiency, state-of-the-art techniques to obtain these functions can be based on manipulations of the symbolic expressions of d z , obtained using the aforementioned recursive formulations. The ideas presented in this paper are more clearly explained making explicit the linearity of d z with respect to the dynamic parameters [26] d This is a standard form of the IDM used in the literature on dynamic parameter estimation, where K zϕ is the singlesample observation matrix. The goal of model simplification by parameter reduction can be defined in this context as eliminating some parameters from ϕ and the corresponding columns of K zϕ so that vector τ z can be approximated as a linear combination of a subset of its columns as where ϕ R is the reduced parameter set and K zϕ R the reduced single-sample observation matrix. For IDM implementation purposes, it is preferable to eliminate the parameters directly from the recursively computed symbolic expressions for d z , that is, approximating it as because in general this explicit linear form K zϕ R ϕ R is computationally less efficient. The computational cost of approximating τ z with the reduced model is obviously smaller. As the functions M zz and δ z are also functions of ϕ, the DDM functions M zz and δ z can be reduced using the same idea. For efficiency purposes, as functions share common subexpressions, they are usually evaluated in a single function call. To illustrate this idea this function will be represented below as [M zz |δ z ](z,ż, ϕ). So, for both IDM and DDM the computational complexity is reduced by symbolically simplifying the terms that are multiplied by the removed parameters in the functions d z (z,ż,z, ϕ) and [M zz |δ z ](z,ż, ϕ), respectively.

Base Parameter Reduction
In the context of dynamic parameter identification of robotic systems, the model is frequently re-parametrized in terms of a smaller independent set of dynamic parameters. One such parametrization is known as the base parameter set [27], or minimum parameter set. This parametrization, which leads to an exact model, can be regarded as a model reduction technique.
The parameter estimation problem is set up by ensuring that the dynamic equations shown in Eq. (5) are satified for a set of observations or estimation set, E = ( i z, iż , iz , i τ z )| i = 1, . . . , n E : where W(E) is the so-called observation matrix for the data set E.
In general, no matter how "exciting" the estimation set is, W(E) is rank deficient, meaning that there are linear dependencies between its columns and, therefore, the set of parameters ϕ is not independent. W(E) can be reordered as where matrix W ϕ R is made up of a chosen 1 set of "independent" (or reduced) columns and matrix W ϕ E is formed by the remaining "dependent" (or excluded) columns. Parameters ϕ can be rearranged accordingly as where ϕ R are the set of "independent" (or reduced) parameters associated with W ϕ R and ϕ E are the set of "dependent" (or excluded) parameters associated with W ϕ E . The columns of W ϕ E can be expressed as linear combinations of the columns of matrix W ϕ R as: In this expression, the columns of matrix β ϕ R contain the coefficients of said linear combinations. Using this relation, Eq. (8) can be rewritten as leading to the base parameter identification problem where is the so-called base or minimum parameter set. In this context, the number of parameters in ϕ R is n ϕ R = rank(W(E)). Note that this parametrization is dependent on ϕ R , the set of independent parameters chosen. Also, for sufficiently exciting data sets (i.e. if E maximizes rank(W(E)) and W ϕ R is well conditioned) ϕ ′ R and β ϕ R become independent of E. In particular β ϕ R becomes dependent only on the geometric parameters of the system. Accordingly, K zϕ can be split as [K zϕ R , K zϕ E ]. Then the dynamic model can be expressed in terms of the base parameters ϕ ′ R , as follows: This equation can be regarded as a parameter-reduced model that reproduces the original one exactly. In [12] this reduced parametrization is applied to the IDM of some serial manipulators to reduce their operation count without compromising their accuracy.

Model parameter reduction
The results in the context of parameter estimation of multibody systems show that frequently not all the base or minimum parameters can be estimated with "sufficient" accuracy. This may be related to the fact that for a given estimation data set, E, the columns of matrix W(E) can be expressed with "sufficient" precision as linear combinations of a subset W ϕ R of n ϕ R , n ϕ R < rank(W(E)), independent columns of W(E). It can be said that the characteristic data set does not sufficiently excite certain system dynamics.
Let ϕ R be the parameters associated with the chosen independent column subset W ϕ R . The associated reduced model will be referred to as the ϕ R -parameter-reduced model, or ϕ R model for short.
Let W ϕ E be the set of removed columns and ϕ E their associated parameters. Proceeding as in the previous section, the columns of W(E) and the set of parameters ϕ can be reordered as W ϕ E can be approximated by its projection onto the column space of W ϕ R as: where Based on the preceding two equations, W(E)ϕ can be approximated as: leading to the approximate estimation problem: The previous relations show that the least square solution to this estimation problem is: where β ϕ R (E) can be obtained using Eq. (17). Eq. (20) can be seen as a generalization of the base parameter concept when the number of parameters used is smaller than the rank of W(E), n ϕ R < rank(W(E)), that is, a generalization of the concept of base parameters to the context of approximate parameter-reduced models. We use the term "approximate base parameter set" or "approximate minimum parameter set" to distinguish it from the exact base parametrization described in the previous section. This equation can be used to compute the values of parameters ϕ ′ R (E) from the original model values ϕ. Additionally, this equation provides information on the influence that the different parameters have on the dynamics. Note that this information is not contained in the bare numerical values of ϕ ′ R (E). The parameter-reduced model selection problem can be formulated as follows: "Given a dynamic model with a known parameter set determine a minimal set of parameters of the model that approximates the characteristic model data with the desired accuracy." This leads to the questions of defining: 1) the characteristic estimation data set, and 2) the error measure.
1) The characteristic estimation data set is system and application dependent, and should characterize the system's entire dynamic range in the considered application. This set, E = ( i z, iż , iz , i τ z )| i = 1, . . . , n E , can be obtained by sampling data from full model simulations.
For example, in the case of fully-actuated multibody systems, trajectories are parametrized for a set of independent coordinates, and an optimization criterion is used to determine the values of their parameters. A comparison of some of the classic criteria for trajectory optimization can be found in [28]. Trajectory parametrizations based on harmonic series [29,30,31] and on polynomials [32,33] have been proposed in the literature. For underactuated systems, dynamic simulations must be performed for a set of representative working simulations. The validation data set, . . , n V , is obtained using the same procedure.
2) For a given ϕ R parameter-reduced model, the following normalized error criterion or prediction error measure for the IDM is proposed: In this expression W ϕ R (E) is the observation matrix of the ϕ R -parameter-reduced model determined using the data set E. Parameters ϕ ′ R (E) are the numerical values for the generalized base parameters, defined in Eq. 20, for the ϕ R model using the data set E. The weighting matrix Σ 1/2 is defined as: where nom(τ z ) represents a vector of characteristic values for the elements of vector τ z . In general, characteristic values based on the estimation data can be used (e.g. standard deviation), although they may be established by other means (e.g. nominal torque value of an actuator).
It is important to note that for function ǫ τ z (ϕ R , E) the argument ϕ R refers to the reduced set of parameters and not to their numerical values. Obviously, the computation of this function requires the determination of the numerical values ϕ ′ R (E). Indexes to characterize the relevance of the contribution of different dynamic forces (Coriolis, centripetal, etc.) to the IDM have been proposed [34], but these are not well suited to the parameter reduction problem on which this work focuses.
In a more programmatic style, the error criteria can be expressed equivalently as: where col(S) arranges the elements of S into a column, and i d z and i d z (ϕ ′ R (E)) are external forces for the i-th data sample predicted by the full and ϕ R parameter-reduced IDM, respectively. Note that, for short, . ϕ are the actual parameter values and ϕ ′ R (E) are the generalized base parameters defined in Eq. 20. In Eq. (21) and (23), the proposed norm can be defined in different ways. The most interesting norms in this context are probably the Euclidean norm (chosen in the examples in this paper), the 1-norm and the infinity norm.
In this paper, this error measure, ǫ τ z (ϕ R , E), is used as the fitness criterion for the choice of parameter-reduced models and for their validation. For example, if the number of parameters, n ϕ R , of the desired ϕ R model is fixed, the best model can be determined as where card(ϕ R ) is the cardinality of the set of parameters ϕ R . As before, note that the argument ϕ R refers to the reduced set of parameters chosen and not to their numerical values. For the desired error level, ǫ τ z (ϕ R , E), the reduced model is considered acceptable if validation and estimation data give similar error levels Analogously, the normalized error criterion or normalized prediction error for the DDM can be defined as where iz and iz (ϕ ′ R (E)) are the accelerations for the i-th data sample predicted by the full and ϕ R parameterreduced DDM, respectively. Although the model selection performed in this article is not based on this error, to enrich the discussion, values of this error measure will be computed for the parameter-reduced models obtained based on ǫ τ z (ϕ R , E).
The problem with this approach is that the number of possible parameter sets to be tried in order to find the best parameter-reduced model with card(ϕ R ) = n ϕ R is huge: If all parameter sets with an arbitrary number of parameters n ϕ R = 1, . . . , n ϕ are tried, a much larger number of candidate sets, 2 n ϕ , is obtained. In general, this excludes a trial-and-error approach to the parameter reduction problem. Therefore, algorithms that find a compromise between the computational cost of the search and the quality of the reduced model are required. To that end, three candidate model selection heuristics are proposed: QR selection in Section 2.3, Backward Elimination in Section 2.4, and Forward Selection in Section 2.4.

QR decomposition-based heuristic
A prototypical QR decomposition with column pivoting of the observation matrix W(E) takes the following form: where matrix WE is a column permutation of W(E), Q is an orthonormal matrix, and R is an upper triangular matrix with diagonal elements of decreasing magnitude. Let r = n ϕ R ≤ rank(W(E)). WE is split as where W ϕ R contains the first r columns of WE, and R rr is an r × r regular matrix. The reduced model is defined by the parameter set ϕ R associated with the columns of W ϕ R . From the previous equation it follows that Based on Eq. (16) matrix W ϕ E can be approximated as: where Eq. (30) is a very good approximation of W ϕ E for a parameter set of r elements. Note that R ee is chosen by the QR algorithm to have the smallest diagonal elements and, therefore, the error in the approximation of W ϕ E , Q E R ee , is nearly as small as possible. The ideas presented here are partly inspired by the QR decomposition-based base parameter determination algorithm described in [27], although the ideas presented there focus on the case of r = n ϕ R = rank(W(E)). A subset selection algorithm with some resemblance to the one presented here, referred to as Orthogonal Least Squares, is described in [23]. The normalized prediction error ǫ τ z (ϕ R , E) corresponding to the parameters associated with the first r columns of WE can be computed using Eq. (21). Given a maximum acceptable error tolerance, tol, and starting with n ϕ R = rank(W(E)), we decrease n ϕ R while ǫ τ z (ϕ R , E) < tol. Then, the simplest parametrization based on the QR heuristics is that corresponding to the parameters associated with the first n ϕ R rows of WE.
A possible criticism of this algorithm is that the parameter order obtained is only dependent on W(E), but not on χ(E). The results presented in Section 3 will prove that the heuristics proposed in the next two sections perform much better than this one.

Backward Elimination heuristic
The Backward Elimination (BE) algorithm is an iterative algorithm. The algorithm starts by using the full set of parameters for the current parameter-reduced model, ϕ R =ϕ. At each iteration, the parameter with the least significant contribution to the normalized error, ϕ e = arg min is eliminated from the current model parameter set, ϕ R = ϕ R ϕ e , until the normalized error criterion is greater than the chosen tolerance, ǫ τ z (ϕ R , E) > tol. This is shown in pseudocode in Algorithm 1 using a more formal algorithmic style. Basically, the parameters that make a smaller contribution to the modeling error are removed one by one. The algorithm would be optimal if the contributions to the error of the different parameters, ϕ, were independent. The issue is that these contributions are not generally independent, so eliminating a given parameter has an effect on the weight that the remaining parameters have for the resulting reduced model. Therefore, the algorithm is likely to perform well but not optimally.
An important feature of this algorithm is that, in addition to W(E), the vector χ(E) has an effect on the order in which the parameters are removed. This information is not considered by the previously introduced QR heuristic, so this may be a plus.

Forward Selection heuristic
The iterative Forward Selection (FS) algorithm is, in a way, the reverse of the previous one. It starts with an empty parameter set for the current parameter-reduced model, ϕ R =∅. At each iteration, the parameter with the most significant contribution to the normalized error, is added to the current model parameter set, ϕ R = ϕ R ∪ ϕ s , until the required tolerance objective is met, ǫ τ z (ϕ R , E) < tol. This is shown in pseudocode in Algorithm 2 using a more formal algorithmic style.
In essence, the parameters that make a larger contribution to the modeling error are introduced one by one. As before, the algorithm would be optimal if the contributions to the error of the different parameters, ϕ, were independent. The issue is that these contributions are not generally independent, so adding a given parameter has an effect on the weight that the parameters already introduced have for the resulting reduced model. Therefore, the algorithm is likely to perform well but not optimally.
Later in this paper, it will be seen that FS performs slightly better than the BE heuristic: at the beginning of BE, the model can be expressed exactly in terms of the base parameter set. As this set is not unique, there are several possible parameter removals that have no effect whatsoever on the error. Therefore, there are no criteria to differentiate among possible parameter eliminations but, still, the choice has an effect on the performance of the resulting parameterreduced model. The FS algorithm never has to face the selection of noncontributing parameters, as these are never added.

Other Forward and Backward heuristics
The algorithms presented in Sections 2.4 and 2.5 are analogous to the forward selection and backward elimination algorithms typically found in the context of model selection [22,23]. The literature suggests the possibility of using more complex step-wise [23] iterations in the forward/backward selection/elimination heuristics. Typically, one or more search iterations are introduced in addition to the initial one.
For example (adapted from [22]), a possible forward selection algorithm with an additional search iteration would read: let ϕ R be the set of parameters at the start of a given iteration, and let ϕ l be the parameter added to this set in the previous iteration. The algorithm starts with an empty parameter set for the current parameter-reduced model, ϕ R =∅, and an empty variable as the one added in the previous iteration, ϕ l = ∅. At each iteration of the forward selection, as before, the parameter with the most significant contribution to the normalized error, is added to the current model, ϕ R = ϕ R ∪ ϕ s . Then, the parameter added in the previous iteration, ϕ l , is removed, ϕ R = ϕ R ϕ l . Next, the parameter with the most significant contribution to the normalized error is searched anew, and added to the set, ϕ R = ϕ R ∪ ϕ s , until the error tolerance is met, ǫ τ z (ϕ R , E) < tol. This is shown more formally in the pseudocode in Algorithm 3.

Algorithm 3 Forward Selection with two iterations
This algorithm is not going be tested on the application examples analyzed in the next section. It has just been presented here to give some perspective to the ideas presented in this paper.

Examples and Simulations
In this article, high/low mobility of a body in a multibody system is defined as the possibility of the body to exhibit large/small excursions in its relative motion space with respect to other bodies and/or with respect to the inertial frame. It plays an important role in the relevance of different dynamic parameters. For example, the ideas presented in [35] show that for each d.o.f lost between two arbitrary bodies, or between an arbitrary body and the inertial frame, the number of base parameters of the system can be reduced. Thus, if some relative motion between two bodies or between a body and the ground has a small amplitude, the relevance of some parameters in the dynamics might be small. Therefore, in these cases, it is expected that precise parameter-reduced models requiring a relatively small number of parameters will be obtained. This is why we propose to test the different parameter reduction strategies on two different multibody systems with different -high and low-characteristic mobility. The results for the selected high-mobility system are presented in Section 3.1 and the results for the low-mobility system are given in Section 3.2. First, the system description, the modeling approach used and the procedure followed to obtain the characteristic estimation, E, and validation, V, data sets are described in detail. Then, the results of applying the proposed heuristics to obtain parameter-reduced models with an arbitrary number of parameters, n ϕ R = 1, . . . , n ϕ , are presented in a way that is appropriate for the discussion. As stated before, the normalized error for the inverse dynamics model, ǫ τ z (ϕ R , E), is the performance measure used to obtain every parameter-reduced model.

High-mobility multibody example
The six d.o.f. PUMA560 serial manipulator, shown in Fig. 1, is the high-mobility multibody system example (HME) selected. The modified Denavit-Hartenberg [36] kinematic parameters and coordinates and the inertia parameter set used to define the model, taken from reference [37], are shown in Tables 1 and 2, respectively.
The Principle of Virtual Power has been used to obtain the dynamic equations. As joint/relative coordinates are used for the modeling, no constraint equations and no Lagrange multipliers are required. The symbolic library in [38] has been used to obtain the IDM and DDM functions, d z (z,ż,z, ϕ R ) and [M zz |δ z ](z,ż, ϕ R ), respectively. This library strongly minimizes the number of algebraic operations using algorithms compatible with O(n 3 ) recursive formulations, with a performance that is on a par with state-of-the-art formulations. In this way, the reported simplifications are

Estimation trajectories
The condition number of the observation matrix κ(W(E)) 2 has been used as the objective function for optimization of the trajectories. The independent coordinates of the robot, z 1 , . . . , z 6 , have been parametrized using a finite Fourier series, as proposed in [30]. A total of 100 sampling data points evenly covering a full period are extracted from each optimized trajectory. In order to obtain a characteristic sample set that sufficiently excites the system under the defined optimization conditions, the sampled data from 10 different optimized trajectories is used for the estimation. The MATLAB Optimization Toolbox function fmincon has been used to perform the actual optimization. The rotation angles and the angular velocity ranges of the actuators have been limited using linear inequality constraints. For reference, the details of the trajectory optimization setup are summarized in Table 3. To give a graphical idea, one of the optimized trajectories is shown in Fig. 2.

Performance of the parameter-reduced models
The performance of the previously proposed algorithms is shown in Figs. 3, 4, 5, 6 and 7. Fig. 3 represents the normalized prediction error for the parameter-reduced inverse dynamic models ǫ τ z , in terms of the number of parameters n ϕ R , using the three proposed heuristics for both the estimation, E, and validation, V, data sets. Fig. 4 shows the number of operations, n op , required to compute the functions d z (z,ż,z, ϕ R ) and [M zz |δ z ](z,ż,z, ϕ R ), in terms of the number of parameters n ϕ R , for the three proposed heuristics. The number of operations, n op , refers to the total number of non-repeated sums, products and function evaluations required. Fig. 5 represents the normalized prediction error ǫ τ z in terms of the number of operations, n op for the same functions and heuristics as in Fig. 4. Fig. 6 compares the normalized torque for the full model and the parameter-reduced model with n ϕ R = 18 obtained with the Forward Selection heuristic for the validation trajectory. Fig. 7 shows the normalized error for the parameter-reduced DDM, ǫz, in terms of the number of parameters n ϕ R , using the proposed heuristics for the estimation, E, and validation, V, data sets.
Finally, for completeness, the parameter ordering corresponding to the different heuristics analyzed is presented in Table 4. For a given heuristics row, the set of parameters corresponding to the parameter-reduced model with n ϕ R parameters are those related to the column indexes 1, . . . , n ϕ R .   Figure 3: HME: ǫ τz (ϕ R , E) and ǫ τz (ϕ R , V) vs. n ϕ R for the QR, FS and BE heuristics.

Low-mobility multibody example
The particular difficulties of parallel mechanism identification are well known. For example, in [40] the problems associated with an automobile suspension system are analyzed. These are frequently related to the presence of closed loops that limit the mobility of some bodies to a significant extent. In this study, the 6-PUS Hexaglide [41], shown in Fig. 8, is the selected low-mobility multibody system example (LME).
A total of 24 generalized coordinates, q, are used to model the system: the heights of the 6 guides, z = [z 1 , . . . , z 6 ], that are the independent/actuated coordinates, 2 relative rotations of the universal joints for each of the 6 arms, and 6 absolute coordinates for the head (the Cartesian coordinates of the head center, and its Euler angles). For reference purposes, the kinematic parameters of the prototype are shown in Table 5. The closed-loop nature of the system requires the enforcement of 18 constraints (3 for each S bar-head joint), leading to a total of 6 d.o.f. As before, the Principle of Virtual Power has been used to obtain the dynamic equations. The inertial parameters of the bars, referenced to the U joint with the carriages, and those of the Head, referenced to the center of the S joints, are given in Table 6. The required constraint enforcement is performed using Lagrange multipliers λ, leading to the set of dynamic equations where φ q = ∂φ ∂q is the constraint Jacobian. This equation can be rearranged as The generalized coordinate vector is ordered as q = d T , z T T , where z is the independent coordinate set (linear actuator coordinates). Accordingly, the Jacobian columns are ordered as φ q = φ d , φ z . Now, noting that is an orthogonal complement of φ q , φ q R = 0, left-multiplying Eqs. (34) and (35) by R T and substituting where γ = φ qq −φ, in Eq. (34), these equations take the form of Eqs. (3) and (5). This leads to the formulation on independent coordinates of the dynamic functions chosen to describe the algorithms introduced in this paper. The state, (q,q), is determined as a function of the independent coordinates and velocities, (z,ż), making use of the constraint equations at the level of position and velocity. As with the previous example, the symbolic library in [38] has been used to obtain the equations. As the parameter reduction procedure does not affect the kinematic model and, therefore, the computation of R, we propose to show the computational complexity of the parameter-reduced models of this example in terms of the number of operations required to compute d q (q,q,q, ϕ R ) and [M qq |δ q ](q,q,q, ϕ R ). These are the matrices directly exported by the symbolic library.

Estimation trajectories
The trajectories are obtained using the same procedure followed for the previous example. The 6 independent coordinates, z 1 , . . . , z 6 , have been parametrized using finite Fourier series. In this case, the characteristic sample set uses data from 25 different optimized trajectories. For each trajectory 400 sample points are extracted. For reference purposes, the details of the trajectory optimization setup are summarized in Table 7. To give a graphical idea, one of the optimized trajectories is shown in Fig. 9.
Performance of the parameter-reduced models For this example, the performance of the proposed algorithms is illustrated in Figs. 10, 11, 12, 13 and 14. Fig. 10 shows the normalized prediction error ǫ τ z , in terms of the number of parameters n ϕ R of the parameterreduced models, for the estimation, E, and validation, V, data sets based on the proposed QR , FS, and BE heuristics. Fig. 11 plots the number of operations, n op , required to compute d q (q,q,q, ϕ R ) and [M qq |δ q ](q,q, ϕ R ), in terms of n ϕ R . Fig. 12 shows the normalized prediction error ǫ τ z in terms of n op for the same functions and heuristics. Fig. 13 compares the normalized torque for the full model and the parameter-reduced model with n ϕ R = 17 obtained with the FS heuristic for the validation trajectory.
Finally, for completeness, the parameter ordering corresponding to the different heuristics analyzed is presented in Table 8.

Discussion
In this section, the results of the previously presented simulations are discussed. An interpretation of the results is given, which aims to draw conclusions that can be extended to a more general MSD context. First, we analyze the reduced parameter models in terms of their accuracy in the range in which the error is negligible. Then, the same analysis is performed for larger reductions with an appreciable effect on precision. Finally, the computational time savings for the different parameter-reduced models are analyzed.

Parameter reduction without loss of accuracy
For a sufficiently exciting estimation trajectory, E, the rank of the observation matrix, r = rank(W(E)), is maximized. In this case, this rank coincides with the number of base parameters, which is 36 and 64 for the analyzed HME and LME, respectively.
For both examples, there are ridges in the plots for ǫ τ z vs. n ϕ R for all the heuristics: QR, FS and BE. These ridges are clearly appreciable in Figs. 3 (HME) and 10 (LME). For the QR heuristic, as expected, the right of this ridge coincides exactly with the number of base parameters. This confirms that the estimation trajectories used in this work are sufficiently exciting and, therefore, appropriate to support the analyses. The models to the right of these ridges will be said to have "full precision", as the small noise appreciable in this zone is related to numerical round-off errors and not to the reduced parameter number.
In the same figures, both for the HME and the LME, it is important to notice that some of the ridges for the FS and BE heuristics are displaced to the left relative to those of the QR heuristic. This can be explained by the fact that the FS and BE heuristics optimize the prediction error ǫ τ z , while QR minimizes the estimation error for the model parameters. For the HME only the FS heuristic achieves full precision for a model with a number of parameters, n ϕ R , that is smaller than r, n ϕ R = r − 1 = 35. In the case of the LME, a reduction in the minimum number of parameters required to achieve full precision is appreciated for both the FS and BE heuristics. The FS heuristic clearly performs best followed by BE, requiring a minimum of n ϕ R = r − 5 = 59 and n ϕ R = r − 2 = 62 parameters, respectively. The fact that the distance of the FS and LS ridges to the left of the QR ridge is larger in the LME than in the HME is a direct consequence of its "low mobility".
These ridges also appear in the plots for ǫz vs. n ϕ R shown in Figs. 7 (HME) and 14 (LME). For all the heuristics, QR, FS and BE, these ridges are placed exactly at the same number of parameters, n ϕ R , as in the ǫ τ z vs. n ϕ R plots discussed in the previous paragraphs. So the same conclusions about the attainable reduction in the number of parameters, while maintaining full precision, are applicable not only to the IDM but also to the DDM.
In the light of the aforementioned similarities, noting that the results in this paper are based on the minimization of ǫ τ z , it seems that a hypothetical model selection based on the minimization of ǫz, instead of the minimization of ǫ τ z , will produce similar results for parameter-reduced models with full precision.

Parameter reduction with small loss of accuracy
Now, attention is focused on the parameter-reduced models on the left of the ridges described above. In this zone, a loss of accuracy which is not occluded by numerical round-off errors can be observed.
Going back to the plots for ǫ τ z vs. n ϕ R shown in Figs. 3 (HME) and 10 (LME), it can be seen that the proposed FS and BE heuristics outperform the QR heuristic, showing errors of a smaller order of magnitude for the same number of reduced parameters. In fact, the QR heuristic cannot be claimed to be a reasonable model selection or parameter reduction method on the left of its ridge. As before, this is explained because the QR heuristic minimizes the error on the estimation of the parameters and not the prediction error ǫ τ z . This puts the relevance of the performance indexes chosen to obtain parameter-reduced models into perspective.
These figures clearly show that the FS and BE heuristics can be used to obtain very significant reductions in the number of parameters while still keeping a respectable level of accuracy. Based on the FS heuristic, errors of about ǫ τ z ≈ 10 −2 can be obtained for the HME and the LME with as few as n ϕ R = 18 and n ϕ R = 17 parameters, respectively. These are the parameter-reduced models selected for discussion throughout this section. See Table 9, which summarizes the relevant data used in the discussion of these examples.
Even larger reductions in the number of parameters can be obtained if precision requirements are relaxed a bit more. For a given n ϕ R , at the left of the corresponding ridge, the FS and FE heuristics give a similar ǫ ϕ R error. Nevertheless, FS consistently performs better than BE in both examples.
In the same zone, the plots for ǫz vs. n ϕ R shown in Figs. 7 (HME) and 14 (LME) depict the same tendency: parameter-reduced models with the same number of parameters as before, n ϕ R = 18 for the HME and 17 for the LME, show acceptable ǫz error values (see Table 9). This error seems to increase rapidly if the number of parameters is further reduced. In both examples, even for the best-performing FS heuristic, errors > 100% for a parameter count of ≈ 10 can be observed. This sharp increase is related to the degradation of the condition number of the mass matrix M zz , whose inverse is required to computez. Note that this error goes up to +∞ when M zz becomes singular. This suggests that in order to improve the performance of the parameter-reduced DDM, ǫ τ z can be substituted as the optimization criterion by ǫz. This may be of particular interest when dealing with parameter-reduced models with a small number of parameters. An interesting benefit of using this error as the optimization criterion is that it does not become singular for underactuated multibody systems.
Again it is seen that, for an equal number of parameters, the FS heuristic shows smaller errors than the BE heuristic. Let us use the LME to explain this difference: the FS heuristic algorithm begins with the full set of parameters, ϕ R = ϕ. When it starts removing parameters, it faces the removal of parameters that have a negligible effect on the dynamics. With n ϕ R between 70 and r = rank(W(E)) = 64, at every reduction step the removal of parameters has no effect at all on the dynamics. There is a large number of parameters that can be chosen to be removed, and the algorithm randomly decides which one because, in this range, ǫ τ z is just numerical noise. But this initial choice of removed parameters affects the parameter-reduced models obtained afterward. On the other hand, the FS algorithm starts with an empty set of parameters. At each step, it always adds parameters with a non-negligible contribution to the dynamics and, even if it is not optimal, it will always capture the leftmost ridge position. As the results show, this makes a difference. A similar reasoning can be applied to the HME. These results are in accordance with results observed for the FS and BE heuristics in other scientific contexts [23].

Operation count savings vs precision
The number of operations in relation to the number of parameters, n op vs. n ϕ R , for the different heuristics is shown in Figs. 4 (HME) and 11 (LME). It is important to note that the criteria used for model selection do not include minimization of the operation count of the reduced models. Nevertheless, an important decay in the number of operations is observed as the number of parameters in the model is reduced. This decay exhibits an almost piece-wise linear profile, with slopes that increase with the number of parameters removed. The figures suggest that there are some advantages to the FS and BE heuristics in comparison with the QR heuristic. It is important to remember that the operation count is dependent on different aspects, such as the parametrization used or the formalism and the way in which it is implemented. As mentioned in Section 3, in both examples state-of-the-art general purpose symbolic formulations [10,38] have been used to obtain the model functions and, therefore, the operation count. In this way, the results presented give a correct perspective on the real savings attainable using the proposed methods.
The normalized prediction error as a function of the number of operations, ǫ τ z vs. n op , for the different heuristics is presented in Figs. 5 (HME) and 12 (LME). Note that the information given in these plots is just a convenient rearrangement of the information given in Figs. 3 and 4 (HME) and in Figs. 10 and 11 (LME), respectively. This is a more useful representation of the data if a balance between the precision and computational requirements of the chosen parameter-reduced model is sought. From these figures it can be observed that, for the previously selected models (FS with n ϕ R = 18 for HME and FS with n ϕ R = 17 for the LME, with a prediction error of ǫ τ z ≈ 10 −2 ), the reductions in the computational cost for computation of the IDM and DDM functions are above 46% in the case of the HME and above 67% in the case of the LME (see Table 9). Table 9: Results for selected HME & LME models: heuristic H, parameter number n ϕ , base parameter number r, reduced parameter number n ϕ R , and full and reduced op. number n op (ϕ) and n op (ϕ R ), n op (ϕ R )/n op (ϕ) "red." and normalized errors ǫ τz (ϕ R ) and ǫz(ϕ R In a more direct way, Figs. 6 (HME) and 13 (LME) compare the normalized values of the actuator moments and forces along the validation trajectory for the full and selected parameter-reduced models (FS with n ϕ R = 18 for HME and FS with n ϕ R = 17 for LME), respectively. As can be seen, the trajectories are barely distinguishable by the naked eye.

Conclusions
In this paper, drawing inspiration from the literature on model selection, parameter reduction methodologies have been proposed and tested in the context of MSD models. These methods require: a) a suitable selection heuristic, b) a characteristic data set representative of the dynamics of the system, and c) a performance measure for the parameterreduced models.
Three different model selection heuristics -QR decomposition-based (QR), Forward Selection (FS) and Backward Elimination (BE)-have been proposed and tested in two multibody systems with different -high (HME) and low(LME)-characteristic mobility. A normalized error for the inverse dynamics model is defined and used as the performance measure in these examples. A detailed description of the systems, procedure and results is given.
The results show that, based on the FS heuristics and the proposed performance measure, an impressive reduction in the number of parameters required is obtained: for a normalized error for the inverse dynamics model of ≈ 1% only 18 parameters out of 49 for the HME and 17 out of 70 for the LME, are required. Corresponding to these parameterreduced models, computational complexity savings for the evaluation of the functions required for the direct and inverse dynamics models are above 36% for the HME and above 67% for the LME. A defined normalized error for the direct dynamics model shows errors of ≈ 5% for these same models.
A detailed description and interpretation of the results is given. In general, the FS heuristic, closely followed by BE, appears to be the best-performing one for any number of parameters. QR, a good heuristic to reduce the error in the parameter estimation context, is discarded due to inadmissible levels for the proposed normalized errors. The degradation of the direct dynamics seen for small parameter-reduced models is related to the worsening of the condition number of the mass matrix.
For any reduction level, the better performance of the methods for the LME is seen to be related to its low mobility. Applications of MSD for low-mobility systems are likely to benefit especially from the methods presented. Important cases in this context are railway and automotive dynamics. For these, the methods proposed can offer very large reductions in the number of parameters, as they present low mobility. This leads to more robust parameter estimation procedures and much smaller computational costs of the IDM and DDM, offering new possibilities to applications of MSD that require real-time performance.
As an additional conclusion, it has been shown that it is possible to obtain analytical expressions that determine the value of the estimated reduced set of parameters as a function of the parameter values of the original model. These analytical expressions can be considered a generalization of the base parameter concept to the case of approximate parameter-reduced models.

Future work
Two lines have been identified as having research potential: first, using the normalized error for the direct dynamics as the optimization criterion for further improving the quality of direct dynamics parameter-reduced models with a small number of parameters. Second, within computational resources, widening the search space of FS-style heuristics to obtain even better results.
With regard to the second line, it is interesting to mention the important case of multibody systems models exhibiting symmetries. The space search can be reduced enormously if characteristic data sets exhibiting the same symmetry are used and if dynamic parameters with a symmetric dynamic role are selected or eliminated in groups. This will generally make it possible to obtain better results, or even optimal results in some cases.