Kinematic and inverse dynamic analysis using mixed and fully Cartesian coordinates with a generic rigid body

Research paper Kinematic and inverse dynamic analysis using mixed and fully Cartesian coordinates with a generic rigid body S´ ergio B. Gonçalves a , Ivo Roupa b , Paulo Flores c , Miguel Tavares da Silva a,* a IDMEC, Instituto Superior T´ ecnico, Universidade de Lisboa, Lisboa, Portugal b ITI/LARSyS, Instituto Superior T´ ecnico, Universidade de Lisboa, Lisboa, Portugal c CMEMS‑UMinho, Department of Mechanical Engineering, University of Minho, Guimar˜ aes, Portugal ARTICLE INFO Keywords: Fully Cartesian coordinates Generic rigid body Mixed coordinates Kinematic analysis Inverse dynamics Advanced education ABSTRACT The propagation of errors along a kinematic chain, caused by using drivers computed from noisy data, can affect the accuracy of the kinematic and dynamic outcomes. Minimizing the effect of such errors is crucial, particularly when traditional smoothing techniques prove to be ineffective. This work expands the multibody formulation with Fully Cartesian Coordinates and a Generic Rigid Body (FCC-GRB) to the inverse dynamic analysis of spatial mechanical systems using mixed coordinates (MC). This method considers the incorporation of angular variables, enabling the determination of the kinematic consistent positions that best fit the reference data, while simultaneously computing the joint angular drivers. The accuracy and computational performance of the formulation are evaluated using both numericaland optimization-based methods in the study of two mechanisms guided with perturbed data. The results show that implementing an MC methodology with FCC-GRB can be easily performed without compromising the theoretical foundations of the classical formulation. This approach efficiently computes both the positions and drivers of the model simultaneously, avoiding the propagation of errors along the kinematic chain. Numerical methods based on the Newton-Raphson algorithm generated positions closely matching the reference data, while optimization-based methods ensured a stricter fulfillment of the topological constraints. 1. Introduction The inverse dynamic analysis of multibody mechanical systems typically requires prior computation of kinematic positions from noisy or missing experimental data [1,2]. A classic example is the biomechanical analysis of human motion. Typically, when acquiring experimental data, a set of reflective markers is placed on the subject’s skin at specific locations, often corresponding to bony eminences [3–5]. Due to the relative motion of soft tissues, such as the skin, muscles, tendons, and ligaments, in relation to the bones, a phenomenon usually referred to as soft tissue artifact (STA), the acquired data is prone to experimental errors. These can affect the quality of the measured outcomes, making it challenging to use the data at a clinical level [6–10]. This artifact is particularly difficult to remove from experimental data without relying on specialized real-time medical imaging techniques (e.g., [11–13]), as it often has the same frequency as the movement itself and is task-dependent. For instance, its magnitude, frequency, and timing may vary * Corresponding author. E-mail addresses: [email protected] (S.B. Gonçalves), [email protected] (I. Roupa), [email protected] (P. Flores), [email protected] (M.T. Silva). Contents lists available at ScienceDirect Mechanism and Machine Theory journal homepage: www.elsevier.com/locate/mechmt https://doi.org/10.1016/j.mechmachtheory.2025.106080 Received 24 March 2025; Received in revised form 8 May 2025; Accepted 15 May 2025 Mechanism and Machine Theory 214 (2025) 106080 Available online 5 June 2025 0094-114X/© 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). Nomenclature Head Abbreviations 3D Three-Dimensional AR Angular Relation CNO Constrained Nonlinear Optimization CoM Center of Mass DoF Degree of Freedom ED Euclidean Distance EoM Equations of Motion FCC-GRB Fully Cartesian Coordinates with Generic Rigid Body KA Kinematic Analysis LS Least Squares MCR Mixed Coordinates Relation MO Multi-Objective NRM Newton-Raphson Method STA Soft Tissue Artifact TC Tracking Constraint TOC Tracking Orientation Constraint WLS Weighted Least Squares Symbol (Latin) 03, I3Null and identity matrices CP iConstant transformation matrix of a generic point P with respect to body i Cs iConstant transformation matrix of a generic vector s with respect to body i Cs1s2Constant transformation matrix that relates the generic vectors s1 and s2 F, FObjective function of the optimization problem FqJacobian matrix of the objective function F F∗Optimization goals of the multi-objective goal attainment problem gGeneralized forces of the system gΦGeneralized internal forces MMass matrix of the system nb, nc, nd, nTC Number of bodies, kinematic constraints, generalized angular variables and tracking constraints of the system PGeneric point P P∗Reference point P∗ OiOrigin of the local reference frame of body i q, ˙ q, ¨ qGeneralized positions, velocities and accelerations of the system q3vi, ˙ q3vi, ¨ q3viPosition, velocity and acceleration of the fully-defined equivalent vector of the reduced rigid body i qlb, qub Lower and upper bounds of the generalized coordinates for the optimization qnSolution a of the NRM for iteration n rP,˙ rP,¨ rPPosition, velocity and acceleration vectors of the generic point P in the global reference frame rP∗, ˙ rP∗, ¨ rP∗Position, velocity and acceleration vectors of a reference point P∗in the global reference frame s, ˙ s, ¨ sPosition, velocity and acceleration vectors of a generic unit vector s s∗, ˙ s∗, ¨ s∗Position, velocity and acceleration vectors of a reference unit vector s∗ SConstant transformation matrix that converts the reduced rigid body i in its fully-defined equivalent form tTime ui,˙ ui,¨ ui ,, vi,˙ vi,¨ vi ,wi, ˙ wi, ¨ wiPosition, velocity and acceleration of the generic rigid body vectors u, v and w of body i in the global reference frame Vi,Vi,˙ ViTransformation matrices that convert the reduced rigid body i into its fully-defined equivalent form wnWeight coefficient of the n-th tracking constraint wMO Weight coefficients vector for the MO optimization WNRM-WLS and optimization weight matrix Symbol (Greek) γRight-hand-side vector of the acceleration equations of the system γkRight-hand-side vector of the acceleration equations for the kinematic constraint of type k Y Slack variable of the multi-objective goal attainment problem ε Newton-Raphson iteration tolerance θVector of the generalized angular variables θkk-th angular generalized coordinate of the model S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 2 depending on the movement being analyzed, the individual’s anthropometric characteristics, the body segment and anatomical joint involved, and the muscle activation pattern. As a result, automatic detection becomes difficult, and standard filtering techniques are often ineffective for its attenuation [6,14]. Obtaining high-quality kinematic and kinetic inputs is crucial, as they have a significant influence on the quality of the dynamic outcomes [15]. A common procedure for ensuring the quality of kinematic inputs involves computing the kinematic consistent positions for a given model. Essentially, this procedure determines the positions, velocities, and accelerations of the system that comply with a set of topological constraints, which define mathematically the model, and driving constraints that guide its movement [1]. Methodologies based on multibody system formulations present themselves as an efficient and accurate approach to performing the kinematic and dynamic analysis of complex mechanical models [16]. Two classical approaches are commonly employed to compute the kinematic consistent positions [17]. The first addresses the problem through a forward kinematics approach. The joint angular drivers, which usually describe the angular degrees-of-freedom (DoFs) of the model, are computed from the experimental data in a pre-analysis phase. These drivers are then applied to the model during the kinematic analysis (KA) to determine the position and orientation of all the segments [1,17]. The second method utilizes reference points and vectors to determine the kinematic consistent positions that minimize the differences from these references, following an inverse kinematics approach [17,18]. In order to address the inaccuracies in the experimental data, such as STA, several methods have been proposed. Some of them use specific algorithms designed to mitigate these effects (e.g., [19–21]), while others rely on numerical or optimization methods together with multibody-based models (e.g., [7,9,10,22,23]). The choice of the most appropriate method depends on the desired level of accuracy, movement in analysis, or computational performance requirements. Several multibody formulations are available in the literature, presenting different advantages and target applications [17,24–31]. Recently,based on the work of Gameiro et al. [32], Roupa et al. [28] and Gonçalves et al. [29] presented the theoretical basis of the Fully Cartesian Coordinates with a Generic Rigid Body (FCC-GRB) formulation for the planar and spatial analysis of multibody systems. Due to the type of coordinates used to describe the system and the kinematic structure adopted to model the segments, this formulation preserves some of the main characteristics reported for the most common global formulations, namely the natural coordinates [17,27, 33] and Cartesian coordinates [24–26]. Similarly to the natural coordinates formulation, in the FCC-GRB formulation, the multibody system is described using only rectangular coordinates of points and unit vectors, implying that no angular-related variables are required to describe the orientation of the bodies. Furthermore, the kinematics of any point or vector belonging to the model can be expressed as a linear transformation of a set of constant transformation matrices (matrices C and V) and the generalized coordinates of the system. These features imply that the constraint equations for the most common kinematic constraints present, at the most, a quadratic dependency on the generalized coordinates of the system, and, consequently, the respective contributions to the Jacobian matrix and right-hand side vectors of velocity and acceleration are quadratic, linear, constant, or null [28,29]. Conversely to the natural coordinates formulation, in which the definition of the bodies depends on the topology of the system, the rigid bodies are defined with a predetermined kinematic structure composed of one reference point and two or three unit vectors. Moreover, if the reference point is located at the center of mass (CoM) of the body, the rigid body mass matrices become diagonal with their entries equal to the inertial properties of the segments being modeled [28,29]. The use of a generic rigid body approximates the proposed formulation to the Cartesian coordinates approach. Hence, some of the advantages often associated with this formulation, namely the easier systematization of the modeling procedure, and the straightforwardness and high physical meaning of the mass matrices, also apply to the FCC-GRB formulation, simplifying its computational implementation [28,29]. Some authors explored a hybrid multibody formulation based on natural coordinates, which also uses angular variables [34,35]. This approach, originally proposed by Jal´ on and Bayo [17], considers that the set of generalized coordinates includes not only the coordinates that define the model topology but also the angular drivers that guide the system, allowing for the simultaneous determination of the consistent positions and orientations of the model and their respective angular DoFs. Due to the hybrid nature of the coordinates that compose the system, this formulation is commonly referred to as mixed coordinates (MC) [17,34,35]. The modeling approach with MC is particularly useful when performing the kinematic analysis of systems that require the prior calculation of input drivers from experimental data, which may be subject to measurement errors. By minimizing the distance between the experimental data and the corresponding elements of the mechanical model (e.g., points, vectors, orientations), this methodology θ∗ kReference angle for the k-th angular generalized coordinate of the model θs1s2Angle between two generic unit vectors s1 and s2 θ∗ s1s2,˙ θ∗ s1s2,¨ θ∗ s1s2Angular position, velocity and acceleration of a reference angle θ∗between two generic unit vectors s1 and s2 λLagrange multipliers of the system λkLagrange multipliers for the kinematic constraint of type k ν Right-hand-side vector of the velocity equations of the system ν kRight-hand-side vector of the velocity equations for the kinematic constraint of type k τ int kInternal joint moment of force associated to joint k Φ, ˙ Φ, ¨ ΦVector of the kinematic constraints and respective first derivative and second derivative with respect to time ΦqJacobian matrix of the system ΦkKinematic constraint equation of type k Φk qContributions to the Jacobian matrix from the kinematic constraint of type k S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 3 enables to find the kinematic consistent positions that best fit the input data over time. Furthermore, as the entire position is minimized at the same time, small errors due to data inconsistency are not propagated along the kinematic chain, as in the case of the kinematic analysis with drivers. These advantages are particularly suitable for the biomechanics field, where experimental data is often contaminated with STA, which can result in marker displacements of a few centimeters from their original position [8,9,36]. The use of an MC approach allows for attenuation of the effect of these artifacts, while simultaneously computing the joint angular drivers directly from the experimental markers. In addition, this step eliminates the need for pre-processing procedures that involve the calculation of the three-dimensional (3D) angles from experimental data, which is often a complex issue, especially for users with less expertise in the topic. Thus, the present work examines the use of an MC approach in tandem with an FCC-GRB formulation to perform the kinematic and inverse dynamic analysis of spatial multibody systems. Three different approaches to solving the kinematic equations are considered, namely the Newton-Raphson iterative method (NRM) together with a weighted least squares (WLS) minimization, a constrained nonlinear optimization (CNO) and a multi-objective (MO) goal attainment optimization, to evaluate their influence on the accuracy of the reconstruction of the kinematic consistent positions from the experimental data, the violation of the kinematic constraints and computational time. With the purpose of assessing its applicability, the FCC-GRB formulation with an MC approach is applied in the analysis of two mechanisms, considering, as inputs, data with and without perturbations. According to Roupa et al. [28] and Gonçalves et al. [29], one of the major advantages of modeling multibody systems with an FCC-GRB formulation is its ease of implementation and modeling of complex mechanical systems. This fact implies that, in addition to the classical applications usually attributed to the global multibody formulations, the FCC-GRB formulation is suitable for advanced education purposes. As the incorporation of the angular variables does not significantly change the modeling procedure, the straightforwardness reported to the planar and spatial formulation is still applicable in the mixed coordinates variation. Consequently, the steps required to incorporate angular coordinates into the formulation are presented in detail, so that, together with the planar and spatial FCC-GRB studies [28,29], this work can serve as an educational tool for supporting the teaching of multibody-related topics. 2. Formulation 2.1. Fundamental aspects of the FCC-GRB formulation In the FCC-GRB formulation, multibody systems are modeled using a set of rigid bodies defined with a pre-determined kinematic structure. This approach follows the methodology used in the Cartesian coordinates formulation, presenting advantages in the systematization of the modeling procedure and definition of the system mass matrices. However, contrary to Cartesian coordinates formulation case, which uses angular-related coordinates, the mathematical structure adopted for defining each body uses only rectangular coordinates of points and vectors. Specifically, the structure adopted for the rigid body considers the use of one vector that defines the reference point (rOi), usually located at its center of mass, and three non-coplanar unit vectors (ui, vi, wi), as schematically represented in Fig. 1a. Therefore, the vector of the generalized coordinates for the generic rigid body i (qi) is defined as [29] qi={rT OiuT ivT iwT i}T(1) Vector rOi allows to describe the translation movements of the rigid body i, while the vectors ui, vi, wi provide its orientation. Since the three unit vectors define a vector basis in space, it is possible to compute the kinematics of any generic point (P) or vector (s) belonging to the body directly from its generalized coordinates using a set of transformation matrices, as presented below [29] Fig. 1. Representation of the kinematic structure of a generic rigid body i with a generic point P and vector s belonging to it, modeled in its: a) fullydefined form; b) reduced form. S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 4 rP=CP iqi(2) and s=Cs iqi(3) where rP and s represent, respectively, the position vector for the generic point P and vector s in the global reference frame, and CP i and Cs i the respective transformation matrices. As it is possible to establish a vector basis using orthogonalization methods with two non-collinear vectors, the definition of the generic rigid body can be reduced to two unit vectors (see Fig. 1b) qi={rT OiuT ivT i}T(4) The reduced definition, given by Eq. (4), has the advantage of requiring less generalized coordinates and kinematic constraints than the fully-defined form, reducing the dimensions of the problem to solve. However, the definition of the kinematic constraints that describe the topology of the system requires the computation of the fully-defined equivalent vector (q3vi) from the reduced form. This procedure is performed using a set of transformation matrices S and V, which presents an explicit dependency on the generalized coordinates of the system, such as [29] q3vi=SViqi(5) with S= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ I3030303 03I30303 0303I303 030303 1 2I3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦(12×12) (6) and Vi=⎡ ⎢ ⎢ ⎣ I30303 03I303 0303I3 03− vi ui ⎤ ⎥ ⎥ ⎦(12×9) (7) This dependency implies that the degree of the kinematic constraints increases, and, consequently, the respective contributions to the Jacobian matrix (Φq). Moreover, the rigid body mass matrix becomes dependent on the generalized coordinates, also generating velocity-dependent inertial forces [29]. Thus, considering a mechanical model composed of nb bodies, the position and orientation of the multibody system is algebraically defined as Fig. 2. Representation of the kinematic structure and the respective DoFs of a revolute joint defined between vector s1 of rigid body i and vector s2 of rigid body j and of a spherical joint defined between vectors s3 and s4 of rigid bodies j and vectors s5 and s6 of rigid body j +1 modeled using an FCC-GRB formulation with MC. S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 5 q={qT 1⋯qT i⋯qT nb }T(8) where q is the vector of generalized coordinates of the entire system, including all the bodies defined in the fully-defined and reduced forms. 2.2. Mixed coordinates in an FCC-GRB formulation When including an MC approach into the FCC-GRB formulation, the angular drivers (θk) that guide the model are also treated as generalized coordinates of the system to be solved (see Fig. 2). This approach has the advantage of allowing for the simultaneous computation of the coordinates and joint angles of the model that best fit the experimental data, without requiring a pre-processing step to calculate the drivers of the model. By minimizing the distance between the elements of the model and the equivalent elements of the experimental data, the method mitigates the effect of the experimental errors in the kinematic reconstruction of the movement, enabling to find the optimal model configuration that, on a given timestep, best-fits the experimental data. Mathematically, the MC approach involves the addition of a set of angular variables (θ) to the vector of the generalized coordinates q={qT 1⋯qT i⋯qT nb θ1⋯θk⋯θnd }T(9) where θk is the value of the angular displacement of the k-th driver of the model. Considering this approach, the number of angular variables to add should be equal to the number of DoFs to model using this approach (nd). It is worth mentioning that the total number of angular DoFs for the model can be higher, being the other DoFs guided using the traditional approach with kinematic constraint drivers as described in section 2.4.1. The inclusion of the angular coordinates increases the total number of generalized coordinates of the system. Therefore, a set of angular-based constraint equations needs to be added to the vector of the kinematic constraints of the system (Φ) to express the topological dependencies between the FCC-GRB generalized coordinates and the angular coordinates (see Section 2.3). Moreover, additional linear tracking constraints need to be incorporated into the vector Φ to fit the model to the experimental data, effectively guiding the model throughout this process. It should be noted that, when analyzing (bio)mechanical systems, the reference data typically represent the position or orientation of relevant points or vectors in the model. Hence, two main types of tracking constraints are commonly used, namely the tracking constraint (TC) for guiding points and the tracking orientation constraint (TOC) for guiding vectors (see Section 2.4.2). 2.3. Kinematic constraints – Topologic definition As it happens in other multibody formulations, the modeling of mechanical systems with FCC-GRB requires the mathematical definition of their topology. Gonçalves et al. presented a detailed description of the kinematic constraint equations needed for a comprehensive representation of a mechanical system in 3D, using the most common kinematic joints [29]. Consequently, the present work focuses only on the kinematic equations required to fully integrate the MC into the FCC-GRB formulation. 2.3.1. Mixed coordinates relation for a fully-defined body representation From the topological point of view, the use of an MC approach requires the explicit definition of the geometric relations between the generalized coordinates of the FCC-GRB formulation and the newly introduced angular coordinates. This procedure can be described by one algebraic constraint equation that imposes a given angle between two unit vectors belonging to each of the bodies related by the DoF. Let one consider two generic unit vectors s1 and s2 belonging respectively to rigid bodies i and j (see Fig. 2), the mixed coordinates relation (MCR) is described in its homogenous form as ΦMCR(qi,qj,θk)=sT 1s2−cosθk=0 (10) where θk is the angular variable that describes the angular displacement between vectors s1 and s2 (θs1s2). By applying the linear transformation expressed in Eq. (3), the previous equation can be described in terms of the generalized coordinates of the system as ΦMCR(qi,qj,θk)=(Cs1 iqi)T(Cs2 jqj)−cosθk=0 ΦMCR =qT iCs1s2qj−cosθk (11) with Cs1s2=[Cs1 i TCs2 j](12 ×12)(12) Equation (11) expresses a quadratic relation between the generalized coordinates of bodies i and j. However, contrary to the premises of the FCC-GRB formulation, it also includes nonlinear terms that are dependent on the generalized angular variables. Therefore, the contributions to the Jacobian matrix (Φq), comprising the partial derivatives of the kinematic constraint equations with respect to the generalized coordinates of the system, include both the linear terms associated with the generalized coordinates of the two constrained bodies and an additional contribution associated with the angular variable θk, such that S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 6 ΦMCR q={qT jCs1s2T ⏞⏟⏟⏞ ΦMCR qi qT iCs1s2 ⏞⏟⏟⏞ ΦMCR qj sinθk ⏞⏟⏟⏞ ΦMCR θk}(1×25) (13) Due to the constant nature of matrices C, matrix Cs1s2 is also constant and time independent [29]. Moreover, the generalized angular variables do not present an explicit dependency on the time, indicating that Eq. (11) represents a scleronomic constraint. Thus, the contributions to the right-hand side vector of the velocities of the system ( ν ) are null ν MCR =0 (14) In turn, the quadratic and transcendental nature of the MCR condition implies that the contributions to the right-hand side vector of the accelerations (γ) include quadratic terms that depend on the generalized velocities of bodies i (˙ qi) and j (˙ qj), as well the generalized angular velocities (˙ θk), such that γMCR = − 2˙ qT iCs1s2˙ qj−˙ θ2 kcosθk= − 2˙ sT 1˙ s2−˙ θ2 kcosθk(15) It should be noted that Eq. (11) represents one kinematic constraint equation that describes one rotational DoF associated with the kinematic joint defined by bodies i and j. Consequently, to fully describe the rotational DoFs associated with a specific kinematic joint, it is necessary to include, at least, one MCR constraint equation per DoF in the vector of the kinematic constraints, and the respective contributions to the Jacobian matrix and right-hand side vectors of velocity and acceleration. Equations (13) to (15) can be simplified when the vectors that describe the MCR condition are directly rigid body vectors (e.g., ui and vj). In this case, the definition of the MCR condition in its homogenous form is given by ΦMCR(qi,qj,θk)=uiTvj−cosθk(16) This simplicity in the evaluation of the MCR kinematic constraint is also translated into the calculation of the respective contributions to Φq, ν and γ ΦMCR q={vT j ⏞⏟⏟⏞ ΦMCR ui uT i ⏞⏟⏟⏞ ΦMCR vj sinθk ⏞⏟⏟⏞ ΦMCR θk}(1×25) (17) ν MCR =0 (18) γMCR = − 2˙ uT i˙ vj−˙ θ2 kcosθk(19) 2.3.2. Mixed coordinates relation for a reduced body representation As previously mentioned, a reduced approach can be adopted to model the rigid bodies. In this case, the kinematic constraint equations needed to include the MC are similar to those used in the fully-defined case. However, the linear transformation presented in Eq. (5) should be applied to calculate the fully-defined equivalent vector, resulting in the following expression ΦMCR(qi,qj,θk)=(Cs1 iSViqi)T(Cs2 jSVjqj)−cosθk=0 ΦMCR =qT iSViCs1s2SVjqj−cosθk (20) The dependence of the matrices V on the generalized coordinates of the system implies an increase in the degree of the constraint equations. Nevertheless, as discussed in [29], the evaluation of the MCR constraint is highly efficient, as it only considers the algebraic multiplication of linear matrices and vectors dependent on the generalized coordinates of the system. This aspect is also valid for the contributions to the Jacobian matrix and right-hand side vector of accelerations ΦMCR q={qT jVT jSTCs1s2TVi ⏞⏟⏟⏞ ΦMCR qi qT iVT iSTCs1s2Vj ⏞⏟⏟⏞ ΦMCR qj sinθk ⏞⏟⏟⏞ ΦMCR θk}(1×25) (21) γMCR = − (qT iVT iSTCs1s2˙ Vj˙ qj+qT jVT jSTCs1s2T˙ Vi˙ qi+2˙ qT iVT iCs1s2Vj˙ qj)−˙ θ2 kcosθk γMCR = − (sT 1Cs2 j˙ Vj˙ qj+sT 2Cs1 i˙ Vi˙ qi+2˙ sT 1˙ s2)−˙ θ2 kcosθk (22) where ˙ Vi is the time derivative of the transformation matrix V for the body i defined following the methodology presented in [29]. Since Eq. (20) maintains its scleronomic nature, the contributions to the right-hand side vector of velocities for the reduced case remain null. 2.4. Kinematic constraints – Drivers definition In addition to defining the topology of the system, the kinematic analysis of a multibody system requires the description of all the S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 7 DoFs associated with the mechanical model under analysis. Two main approaches can be considered to drive these DoFs, namely a methodology based on translational and angular drivers, or a method based on the minimization of a set of tracking constraints. In the first approach, the angular DoFs need to be calculated from the experimental data, being subsequently fed in the system in the form of kinematic constraints of angular driver type that need to be satisfied for each time frame. The presence of artifacts in the experimental data will be reflected in the corresponding drivers, and, consequently, in the kinematic consistent positions computed during the analysis. A detailed description of the most common driver types used in FCC-GRB formulation can be found in [29]. The second method uses a set of tracking constraints and tracking orientation constraints, which express geometric relations between relevant points and vectors of the model and the respective experimental data, being minimized along the analysis. This method allows for an adaptation of the model to the prescribed data, mitigating the effect of possible experimental errors in the computation of the kinematic consistent positions. It is worth noting that other types of tracking constraints can also be incorporated. For example, the angular relation used to guide the model in the first method could be reformulated as a tracking constraint to be minimized. However, when applied to experimental data acquired using marker-based optoelectronic MOCAP systems, this strategy contrasts with the primary objective of automatically determining joint angular drivers that best reproduce the observed point positions, requiring also a pre-kinematic step to calculate the reference driver value. As this approach falls outside the scope of this study, it will not be explored. Nevertheless, in systems where only the positions of certain points and the joint angles are available, the method could be adapted to include such data, making it possible to compute a solution that fits both the point positions and minimizes the differences in angular drivers. 2.4.1. Angular drivers Within the FCC-GRB formulation, the angular DoFs can be described using an AR condition. This kinematic constraint imposes that the relative angle between two vectors of the system is equal to a prescribed angle. Mathematically, the AR condition is similar to the formulation used to describe the MCR condition. However, it considers that the angular relation is made using the prescribed angle (θ∗ k) instead of the generalized angular variables (θk). Thus, let one consider two generic unit vectors s1 and s2 belonging respectively to rigid bodies i and j defined in their fully-defined form, the AR condition can be described in the homogenous form as ΦAR(qi,qj,t)=sT 1s2−cosθ∗ s1s2(t) = (Cs1 iqi)T(Cs2 jqj)−cosθ∗ s1s2(t) = 0 ΦAR =qT iCs1s2qj−cosθ∗ s1s2(t) (23) where θ∗ s1s2 is the prescribed angle that describes the angular DoF between vectors s1 and s2 along the time. The AR condition expresses a quadratic relation between the generalized coordinates of the system, implying that their contributions to the Jacobian matrix are linear and equal to ΦAR q={qT jCs1s2T ⏞⏟⏟⏞ ΦAR qi qT iCs1s2 ⏞⏟⏟⏞ ΦAR qj}(1×24) (24) Contrary to the MCR condition, the second term of Eq. (23) presents an explicit dependency on the time vector. Hence, the AR condition presents a rheonomic nature, meaning that the contributions to the right-hand-side vector of velocities and accelerations also include terms dependent on the time ν AR = − ˙ θ∗ s1s2(t)sinθ∗ s1s2(t)(25) γAR = − 2˙ qT iCs1s2˙ qj−((˙ θ∗ s1s2(t))2cosθ∗ s1s2(t) + ¨ θ∗ s1s2(t)sinθ∗ s1s2(t))(26) Similarly to the MCR condition, Eq. (23) represents a single kinematic constraint that governs one DoF. Therefore, for each DoF of the model that is guided using an angular driver constraint, it is necessary to incorporate, at least, one AR condition into the vector of kinematic constraints. When the vectors to constrain are rigid body vectors, Eqs. (23) to (26) can be simplified, obtaining ΦAR(qi,qj,t)=uiTvj−cosθ∗ uivj(t) = 0 (27) For this case, the contributions to the Jacobian matrix and right-hand side vectors of velocities and accelerations are identical to those expressed by Eqs. (17) to (19), considering only the terms corresponding to the FCC-GRB generalized coordinates plus the rheonomic terms presented in Eqs. (25) and (26). If a reduced rigid body definition is used for the definition of the model segments, the AR condition can be derived from Eq. (23), considering the computation of the fully-defined equivalent vector as presented in Eq. (5), yielding ΦAR(qi,qj,t)=qT iSViCs1s2SVjqj−cosθ∗ uivj(t) = 0 (28) As for the fully-defined case, the definition of the AR condition in the reduced form is similar to the MCR condition (see Eq. (20)). However, the angular generalized coordinates are substituted by the prescribed angle. Therefore, the contributions to Φ, ν and γ are once again identical to those presented in Eqs. (21) to (22), considering only the terms corresponding to the generalized coordinates of S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 8 the system and the rheonomic terms expressed by Eqs. (25) and (26). 2.4.2. Tracking constraints The formulation of the TC condition is similar to the linear translation driver described in [29]. However, instead of enforcing the distance between a model point and a reference point to be zero along the three axes of the global reference frame, it minimizes the distance between them. Let one consider a generic point P belonging to the rigid body i defined in its fully-defined form, the TC condition can be described in its homogenous form as ΦTC(qi,t) = rPi−rP∗(t) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ rPxi−rP∗ x rPyi−rP∗ y rPzi−rP∗ z ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ =CP iqi−rP∗(t) = 0(29) where rP∗represents the prescribed coordinates of the reference point P* in the global reference frame. Eq. (29) transforms into three linear constraint equations that guide three DoFs associated with point P, implying that the contributions to the Jacobian matrix are constant and equal to the corresponding matrix C ΦTC q=⎡ ⎢ ⎣CP i ⏞⏟⏟⏞ ΦTC qi⎤ ⎥ ⎦(3×12)(30) It should be noted that the contributions can be assembled into the Jacobian matrix of the system during the kinematic preprocessing steps, eliminating the need for its update during the analysis and, consequently, reducing its computational effort. In turn, the rheonomic nature of the TC condition means that the contributions to the right-hand side vectors of the velocities and accelerations are not null, presenting explicit terms depending on the velocity (˙ rP∗) and acceleration (¨ rP∗) of the reference point P*, as presented below ν TC =˙ rP∗(t)(31) γTC =¨ rP∗(t)(32) In a similar way to the kinematic constraints, Eq. (29) can be adapted for the reduced definition of a rigid body using the linear transformation presented in Eq. (5). Thus, let one consider a reference point P* belonging to a rigid body i defined in its reduced form, the TC condition and the corresponding contributions to Φ, ν and γ are given by ΦTC(qi,t) = CP iq3vi−rP∗(t) = CP iSViqi−rP∗(t) = 0(33) ΦTC q=⎡ ⎢ ⎣CP iVi ⏞⏟⏟⏞ ΦTC qi⎤ ⎥ ⎦(3×9) (34) vTC =˙ rP∗(t)(35) γTC =¨ rP∗(t) − (CP i˙ Vi˙ qi)(36) The use of the matrix Vi increases the degree of the kinematic constraint, implying that its contributions to the Jacobian matrix become dependent on the system generalized coordinates, and, consequently, it requires its update in each time step. The tracking constraints can also be applied to map the orientation of a set of guiding vectors. In this case, the kinematic constraint equations in their homogeneous form are similar to those presented for a reference point P*, considering the matrix Cs for a generic vector s (see Eq. (3)) instead of matrix CP. Therefore, let one consider a generic unit vector s belonging to rigid body i and the respective guiding vector s∗, the tracking orientation constraint can be given by ΦTOC(qi,t) = si−s∗(t) = Cs iqi−s∗(t) = 0(37) or in the reduced form ΦTOC(qi,t) = Cs iq3vi−s∗(t) = Cs iSViqi−s∗(t) = 0(38) Mathematically, the TOC formulation is similar to the TC condition, which means that the contributions to the Jacobian matrix (ΦTOC q) are similar to those given by Eqs. (30) or (34), however, using the matrix Cs i. In turn, the contributions to the right-hand side vector of velocities ( ν TOC) and accelerations (γTOC) are equal to those expressed by Eqs. (31) and (32) or (35) and (36), considering the substitution of the transformation matrix CP i and the velocity and acceleration vectors of the reference point by the equivalent entities of vector s∗, respectively ˙ s∗and ¨ s∗. S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 9 kinematic analysis with the NRM (see Figs. C.1 e-f in SM C), with the MO with TC presenting the higher values (~10 –5 m) and the MO with ED with the lower ones (10 –14 to 10 –9 m). Regarding the CNO, no differences were observed between the TC or the ED constraint in the objective function, showing differences to the reference data in the order of 10 –8 to 10 –6 m. 4.1.2. Computational efficiency To begin with, it must be said that all four methods converged quickly to the optimal solution, considering as an initial guess the position of the previous time frame. However, the number of iterations per time step and the computational time depended on the magnitude of the perturbation imposed on the reference data and the method used to perform the kinematic analysis (see Table 1). As expected, the NRM method exhibited a lower number of iterations and faster convergence. In the absence of perturbations and considering a tolerance error of 10 –10 , this method converged in approximately 2.6 s (8.6 ms/frame) for the simpler model, taking around 3 % more time than the kinematic analysis with AD. The increase of the iteration tolerance in four orders of magnitude resulted in a decrease of approximately 1 iteration per time frame and 21.8 % of the computational time. The inclusion of the perturbations in the TC drivers resulted in an increase of the computational time of 8.2 % (9.6 ms/frame) for the normal STA case and 51.7 % (13.5 ms/ frame) for the extreme case. The optimization-based methods for the unperturbed case required more time to find a solution than the NRM with WLS, ranging from 350 % (40.0 ms/frame) for the MO with ED objective function to 1928 % (170.0 ms/frame) in the case of the CNO with interiorpoint algorithm. The addition of the perturbations to the reference data resulted in a similar number of iterations per time frame and computational times of the unperturbed case. A significant difference in the computational performance was observed between the CNO with the default algorithm (interior-point method) and the MO attainment problem, with the latter taking approximately 78 % less time for the unperturbed case. The use of the SQP algorithm allowed to reduce significantly the optimization time for the CNO, reaching an average value of 42 ms/frame. An influence of the objective function can also be found in the MO optimization, with the ED-based function, which generates a lower number of objective functions, taking less time than the one based on the TC condition. 4.2. Five-bar linkage with perturbations in all reference points 4.2.1. Computational accuracy The results obtained for the five-bar linkage with perturbations in all points were similar to those observed in the previous model (see Fig. 5 and Fig. C.2 in SM C). The NRM with LS and WLS allowed to achieve a solution that fulfils with the iteration tolerance defined for the analysis. For the unperturbed case, both methods achieved differences to the reference data in the order of 10 –16 to 10 –12 , while the optimization methods yielded values ranging from 10 –8 and 10 –4 . However, an adaptation of the generalized positions of the model points to the noisy reference points was observed when only the LS approach was utilized, leading to violations of the topological constraints (maximum values of ΦTop iof 3.9 ×10 –3 for a 7.5 mm perturbation and 1.4 ×10 –2 for a 30.0 mm perturbation). The relaxation of the TC weights allowed for the minimization of this issue, obtaining values in the order of 10 –5 and 10 –4 for the analyses with an iteration tolerance of 10 –6 and 10 –10 (see Table 2). The optimization methods enabled the computation of a set of kinematic consistent positions that fit the reference data, while maintaining the maximum constraint violation under 10 –10 . Although the different algorithms generated slightly different solutions, none directly followed the noisy data used as input. A significant difference in the kinematic patterns was observed between the results generated using parallel and sequential optimization, particularly for the perturbed cases. The inability to provide an initial guess near the solution for the parallel method resulted in a noisy pattern, as the method converged to similar positions that also fulfilled the topological constraints but did not necessarily represent the global minimum (see Fig. 5b and d). Moreover, some convergence problems were found in some time frames, even for the unperturbed case, leading to the appearance of some noisy artifacts in the results (see Table 2 and Fig. C.2 in SM C). Table 1 Normalized computational time, average number of iterations per frame, mean constraint violations per frame and maximum value of the constraint violation for the topological kinematic constraints. Normalized Time [ms/frame] N. Iterations ‖Φ‖Max ΦTop i eP0 P7.5 P30 P0 P7.5 P30 P0 P7.5 P30 P0 P7.5 P30 KA AD 10 –6 6.6 - - 3.3 - - 2.0 ×10 –7 - - 7.9 ×10 –7 - -   10 –10 8.6 - - 4.2 - - 2.8 ×10 –12 - - 2.7 ×10 –11 - - NRM LS 10 –6 7.0 7.4 9.1 3.0 3.3 4.1 1.1 ×10 –7 3.1 ×10 –7 2.3 ×10 –7 7.4 ×10 –7 2.4 ×10 –3 8.9 ×10 –3   10 –10 8.9 11.1 13.5 4.0 5.0 6.3 9.9 ×10 –13 1.9 ×10 –11 2.5 ×10 –11 3.5 ×10 –12 2.4 ×10 –3 8.9 ×10 –3 NRM WLS 10 –6 6.9 7.5 9.0 3.0 3.3 4.0 1.1 ×10 –7 2.4 ×10 –7 2.4 ×10 –7 7.4 ×10 –7 1.2 ×10 –5 5.7 ×10 –5   10 –10 8.8 9.6 13.5 4.0 4.9 6.5 1.0 ×10 –12 1.8 ×10 –11 2.4 ×10 –11 2.4 ×10 –12 1.2 ×10 –5 5.7 ×10 –5 CNO IP 10 –10 171.4 173.1 179.2 24.0 24.1 24.3 5.0 ×10 –11 4.9 ×10 –11 5.2 ×10 –11 9.8 ×10 –11 9.9 ×10 –11 9.9 ×10 –11 SQP 10 –10 42.3 42.1 43.0 3.3 3.3 3.4 1.1 ×10 –12 1.8 ×10 –12 4.3 ×10 –12 7.7 ×10 –12 3.6 ×10 –11 9.6 ×10 –11 MO TC 10 –10 90.8 51.6 51.8 12.3 6.6 5.8 5.4 ×10 –11 9.1 ×10 –12 7.1 ×10 –12 9.9 ×10 –11 9.2 ×10 –11 9.9 ×10 –11 ED 10 –10 40.0 41.0 43.3 4.9 5.0 5.1 5.4 ×10 –12 6.6 ×10 –12 9.6 ×10 –12 9.9 ×10 –11 9.6 ×10 –11 9.7 ×10 –11 S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 16 Fig. 5. Comparison between the KA results and reference data for P4 of the five-bar linkage model with a maximum perturbation of 52 mm for all points considering the NRM (left) and optimization algorithms (right): Top (a and b) – x coordinate; Middle (c and d) – y coordinate; Bottom (e and f) – z coordinate. S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 17 4.2.2. Computational efficiency The results for the more complex model with perturbations in all reference points show that the NRM is significantly faster than the best optimization case. Contrary to the simpler model, where the computational times were similar, the increase in model complexity resulted in a 57.7 % increase in computational time for the NRM with WLS with an iteration tolerance of 10 –6 compared to the kinematic analysis with AD. This difference was significantly higher when the iteration tolerance was reduced to 10 –10 , leading to an increase of 79.4 % in the average time per frame, even though an equivalent number of iterations per frame was observed. For the unperturbed case, the CNO with SQP was faster than the MO with both TC and ED objective functions. This trend reverted when the perturbations were added to the input data, with the MO with TC function converging faster. This inversion can be understood by the lower number of iterations and function evaluations that this method requires to converge. Moreover, with the increase of the perturbation, the normalized time per frame reduced for the MO cases. The use of a parallelization strategy with six pools resulted in a reduction of the computational cost, ranging from five times faster in the unperturbed case to two times in the maximum perturbation case. The lower efficiency of parallelization for the more perturbed cases can be explained by the reduction of the relative difference in the number of iterations required to find a solution between the sequential and parallelization cases. In contrast, the reduction for the CNO was approximately two times for all cases. In general, using one additional tracking point led to an increase in computational cost. This difference was more visible in the optimization analyses, with increases ranging from 3.0 ×10 –2 % for the MO with TC and lower perturbation to 43.4 % for the MO with ED and higher perturbation. The results for the NRM with WLS and an iteration tolerance of 10 –10 showed increases ranging from 2.5 % to 6.8 %. 4.3. Five-Bar linkage with perturbation on P 4 and P 9 To simulate the effect of STA in markers more prone to this artifact and evaluate its influence on the propagation of errors along the kinematic chain, a third condition, where only P 4 and P 9 were perturbed, was tested, revealing significant differences between methods (see Fig. 6 and Fig. C.3 in SM C). The results suggest that both the NRM-WLS and the optimization methods can correct the perturbations introduced earlier in the kinematic chain. The positions achieved for the last point of the open chain closely matched the reference position. In contrast, the kinematic analysis with AD and perturbed drivers exhibited substantial discrepancies in terms of the kinematic patterns, reflecting the noise introduced in the data. This difference is clearly evident in the Euclidean distance plot, which showed a maximum distance of 0.05 m for the smaller perturbed case and 0.21 m for the maximum perturbed case. The results indicate that the NRM with WLS tends to generate smaller differences compared to the optimization methods; however, it comes at the cost of not ensuring such strict compliance of the kinematic constraints. Fig. 7 presents one of the joint angular drivers computed for P4 using the MC approach. The results show that even when introducing perturbations in the point that defines this joint and the nearest one P 9 , which assists in tracking the segment longitudinal rotations, the method can compute a consistent angle that matches the reference driver, while reducing the noise. For the normal STA values, the angle pattern follows the reference without significant oscillations (see Fig. 7a). On the other hand, for the maximum perturbed case, although the computed angle follows the expected trend, it still exhibits oscillations throughout the analysis, which are more noticeable in the MO-ED case (see Fig. 7b). It is worth mentioning that the magnitude of these oscillations is significantly lower than what is observed when the drivers are computed from the perturbed data. The comparison between the MC approach and the use of smoothed noisy angular drivers is presented in Fig. 8. The results show that the MC method using NRM-WLS leads to position and velocity profiles closer to the reference data than any of the smoothed driver Table 2 Normalized computational time, average number of iterations per frame, mean constraint violations per frame and maximum value of the constraint violation for the topological kinematic constraints. Normalized Time [ms/frame] N. Iterations ‖Φ‖Max ΦTop i eP0 P7.5 P30 P0 P7.5 P30 P0 P7.5 P30 P0 P7.5 P30 KA AD 10 –6 16.2 - - 3.6 - - 8.2 ×10 –8 - - 1.1 ×10 –6 - -   10 –10 18.4 - - 4.0 - - 6.6 ×10 –12 - - 7.8 ×10 –11 - - NRM LS 10 –6 24.8 36.2 36.2 3.7 5.5 5.6 5.7 ×10 –8 3.2 ×10 –7 3.3 ×10 –7 1.1 ×10 –6 3.9 ×10 –3 1.4 ×10 –2   10 –10 26.10 56.7 57.6 4.0 8.6 8.8 1.2 ×10 –11 3.4 ×10 –11 3.4 ×10 –11 4.3 ×10 –11 3.9 ×10 –3 1.4 ×10 –2 NRM WLS (8 ref pts) 10 –6 25.30 38.8 50.8 3.8 5.8 7.8 3.1 ×10 –8 3.4 ×10 –7 4.2 ×10 –7 9.8 ×10 –7 7.0 ×10 –5 2.5 ×10 –4   10 –10 30.10 60.0 83.5 4.6 9.2 12.9 7.1 ×10 –12 3.1 ×10 –11 4.3 ×10 –11 1.1 ×10 –10 7.0 ×10 –5 2.5 ×10 –4 NRM WLS (9 ref pts) 10 –6 26.6 39.7 41.7 3.8 5.8 8.1 3.2 ×10 –8 3.4 ×10 –7 4.2 ×10 –7 9.0 ×10 –7 7.9 ×10 –5 2.8 ×10 –4   10 –10 31.20 61.5 89.2 4.6 9.2 13.4 5.3 ×10 –12 3.5 ×10 –11 4.2 ×10 –11 1.6 ×10 –10 7.9 ×10 –5 2.8 ×10 –4 CNO (8 ref pts) SQP 10 –10 1199.2 1605.2 1384.7 34.1 46.7 48.1 1.2 ×10 –15 1.3 ×10 –14 2.4 ×10 –13 1.6 ×10 –14 3.8 ×10 –13 3.3 ×10 –11 SQP Par 10 –10 548.2 703.1 554.9 77.0 82.1 83.9 6.3 ×10 –7 4.7 ×10 –6 2.3 ×10 –7 1.9 ×10 –4 1.4 ×10 –3 8.8 ×10 –11 CNO (9 ref pts) SQP 10 –10 1420.3 1786.0 1848.6 38.1 45.9 46.9 1.6 ×10 –15 8.3 ×10 –15 1.1 ×10 –14 2.8 ×10 –14 3.5 ×10 –13 1.0 ×10 –12 MO (8 ref pts) TC 10 –10 1668.5 847.1 504.9 57.7 30.6 22.6 4.4 ×10 –11 1.4 ×10 –11 7.1 ×10 –12 9.9 ×10 –11 2.2 ×10 –10 9.6 ×10 –11 ED 10 –10 2120.7 1016.1 616.4 78.3 40.1 30.7 4.4 ×10 –11 2.6 ×10 –11 3.3 ×10 –11 1.0 ×10 –10 1.0 ×10 –10 1.0 ×10 –10 TC Par 10 –10 343.9 278.8 255.8 78.8 66.1 60.2 3.7 ×10 –11 1.3 ×10 –11 1.3 ×10 –11 1.0 ×10 –10 1.0 ×10 –10 9.9 ×10 –11 ED Par 10 –10 424.7 327.5 278.9 117.3 91.5 75.0 4.8 ×10 –11 2.3 ×10 –11 3.0 ×10 –11 9.9 ×10 –11 1.0 ×10 –10 1.0 ×10 –10 MO (9 ref pts) TC 10 –10 1748.5 847.4 683.1 56.7 28.4 23.0 4.5 ×10 –11 1.5 ×10 –11 9.9 ×10 –12 1.0 ×10 –10 1.0 ×10 –10 9.4 ×10 –11 ED 10 –10 2243.3 1150.1 884.1 76.5 40.6 32.1 4.7 ×10 –11 2.7 ×10 –11 3.1 ×10 –11 1.0 ×10 –10 1.0 ×10 –10 1.0 ×10 –10 S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 18 cases (see Fig. 8a-b), as well as to lower Euclidean distances (see Fig. 8c). Although smoothing the experimental data reduces the influence of the perturbed data on the kinematic results, the Euclidean distances remain approximately one order of magnitude higher than those obtained with the MC approach. For the reference cutoff frequency of 6 Hz, commonly adopted in gait analysis, the average reduction in distance was approximately 21.7 % (see Fig. 8c). Applying a stricter cutoff frequency further reduced both the Euclidean distance and the magnitude of the high-frequency oscillations observed in the joint driver signals (−34.9 %), approximating also the velocity curves of the reference data. However, it should be noted that cutoff values around 3 Hz may be overly restrictive, potentially attenuating meaningful motion components of the signal and thereby compromising the accuracy of both kinematic and dynamic outcomes [46]. The results obtained with the 3 Hz filter further emphasize the advantages of the MC approach over angular-based driving methods, particularly in minimizing the propagation of artifacts along the kinematic chain. Despite reducing noise levels in the joint angles to values similar to those achieved by the MC method (see Fig. 8d), the magnitude of the differences in the point positions remained higher (see Fig. 8c). This outcome stems from the absence of a global fitting process in the angular-based driving methods, implying that perturbations introduced in earlier segments of the chain are directly propagated throughout the system. Furthermore, it should be highlighted that fitting the model to the experimental data enhances not only positional accuracy but also the accuracy of velocity and acceleration patterns, ultimately improving the quality of the dynamic outcomes. Fig. 6. Comparison between the KA results and reference data for P6 of the five-bar linkage model with a maximum perturbation of 52 mm in P4 and P9 considering the AD with perturbed drivers, NRM-WLS and MO with ED: a) x coordinate; b) y coordinate; c) z coordinate; d) Euclidean distance between the KA outcome and reference data. S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 19 Fig. 7. Comparison between the joint angular driver computed using MC and reference data for point P4 in the five-bar linkage model with NRMWLS and MO-ED: a) maximum perturbation of 13 mm in P4 and P9; b) maximum perturbation of 52 mm in P4 and P9. Fig. 8. Comparison between the KA results and reference data of the five-bar linkage model with a maximum perturbation of 52 mm in P4 and P9 considering the AD with perturbed drivers with and without filtering, and NRM-WLS: a) x coordinate of point P6; b) x velocity of point P6; c) Euclidean distance between the KA outcome and reference data for point P6; d) Joint angular driver for point P4. S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 20 5. Discussion The present study expands the fully Cartesian formulation with a generic rigid body to encompass the inverse kinematic and dynamic analysis of multibody systems with Mixed Coordinates and evaluates the accuracy and efficiency of different methods. The adoption of an MC formulation enables the computation of the kinematically consistent positions of the model that better fit the reference data, while simultaneously calculates the angular drivers of the model, which describe the DoFs associated with the kinematic joints. This modeling approach enables the evaluation of the major kinematic outcomes in a single analysis, eliminating the need for additional steps to compute the drivers of the model. Besides simplifying the modeling procedure and preand post-kinematic steps, the presented approach also offers advantages related to the accuracy of the solution. By minimizing the distance between the tracking markers and the corresponding points of the model, the method determines a kinematically consistent position that accurately represents the overall position of the model in relation to the experimental data. This methodology addresses the issues commonly associated with the use of angular drivers, namely the propagation of errors along the kinematic chain due to the computation of drivers from noisy or inconsistent data. Moreover, the use of different weights can force the model to follow points that are less prone to experimental errors or are more relevant to the analysis. It should be noted that this study specifically focused on open-chain models, where error propagation is typically more pronounced due to the absence of topological constraints to enforce positional consistency along the chain. Nonetheless, since the MC method with FCC-GRB maintains the same principles of the global FCC-GRB formulation, it does not face the limitations commonly associated with recursive formulations in the analysis of closed-loop systems [29,47], and, consequently, it can be applied to analyze such systems without requiring significant modifications. From the implementation point of view, incorporating the angular coordinates into the formulation is a straightforward process. The angular DoFs of the system, which are traditionally modeled using angular drivers, are here treated as angular variables that are evaluated during the process of computing the kinematically consistent positions. The AR kinematic constraints, which describe the angular relation between two vectors of the model, are replaced by an MCR condition. In fact, the two conditions are similar, being the rheonomic term of the AR constraint equations substituted by a similar term dependent on the generalized coordinates of the system. This relation implies that the inverse dynamic analysis of a system defined using MC is also a straightforward process. The simplicity of this procedure is related to the fact that one system modeled with MC and FCC-GRB can be easily converted to a classical FCC-GRB model, by transforming the angular generalized variables into angular drivers and by removing the tracking constraints. Therefore, the physical meaning of the internal forces computed using the Lagrange method is the same, not requiring additional steps besides the ones already performed in the classical formulation. It is worth mentioning that although this work details the implementation of a spatial MC formulation based on FCC-GRB, this approach can be applied in other global formulations, considering the specific characteristics of each formulation [17,34,35]. Moreover, transitioning to a 2D FCC-GRB formulation can be easily performed by applying the same concepts as those employed in the spatial formulation, but using the equivalent planar expressions [28]. In this particular case, the key difference lies in the lower minimum number of tracking constraint equations needed to fully describe the system kinematics, as the rigid bodies and kinematic joints have fewer DoFs. The results obtained for the two analyzed models support the previously mentioned advantages. Both the NRM and optimization methods were able to find a kinematically consistent solution, while simultaneously determining the angular drivers of the model. The NRM method was significantly faster than the different optimization conditions tested in this work. This result is a direct consequence of the nature of the NRM algorithm, which exhibits a quadratic convergence near the solution [39,40]. This fact implies that if an initial guess near the solution is given to the iteration algorithm, it will converge quickly to the solution. The use of the NRM method to solve kinematic equations is not without limitations. The classical problems typically associated with the NRM algorithm can occur, including convergence issues when initiating with a poor initial guess or one far from the solution, overshoot effects due to the bad behavior of the function near the solution, or null-derivative functions [39,40,48]. Regarding the optimization methods, the MO attainment goal proved to be the faster solution for the more complex model. The results indicate an influence of the type of objective function utilized. Despite generating more equations than the ED condition, the use of the TC-based objective function resulted in lower computational times for the more complex model, a difference in part explained by the lower number of iterations required to find a solution. This trend was not observed in the simpler model, indicating that the relative differences in the computational times may be strongly affected by the complexity of the model under analysis and the optimization inputs. In fact, the simulations for the cases with higher levels of perturbation converged faster than the ones with lower levels, requiring also less iterations to find an optimal solution. In contrast, the CNO with SQP converged significantly faster than with the IP algorithm. Conversely to the MO case, the increase in the perturbation levels resulted in higher computational times and a higher number of iterations to find a solution. It should be noted that the gradients of the equality constraint equations and objective functions were provided to the optimizer. The use of finite differences significantly increased the simulation time; however, this difference was not quantified in this analysis. In general, the incorporation of the angular coordinates led to an increase in the computational effort for all conditions when compared to the classical FCC-GRB formulation with angular drivers. Nevertheless, it is worth noting that the preand post-processing steps required for computing the drivers are not included in the analysis, and no fitting to the reference data is performed. In contrast, in the MC analyses, both the computation of the generalized coordinates and the angular drivers are carried out simultaneously and are accounted in the computational costs. In terms of accuracy, the kinematic analyses with the NRM method did not yield results as exact as those obtained with the optimization methods. Since the TC and TOC kinematic equations are treated as equality constraints that must be fulfilled, the method generates kinematic positions that, in order to adapt the model to the noisy reference points, may violate some of the constraints. The S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 21 use of different weights for the topological and tracking constraints enables the mitigation of this difficulty. By relaxing the weights associated with the tracking constraints, the method assigns less significance to these equations, prioritizing the fulfillment of the topological constraints. This idea is supported by the results obtained for the analyses with perturbations. The maximum value of the topological constraints violation during the analysis decreased by two orders of magnitude when a weight of 10 –3 was used. This difference can be even increased if lower weights are assigned to the TC equations, or higher weights are given to the topological constraints. This better fitting to the experimental points is also supported by the lower ED to the reference points. Overall, the NRM method presented lower values than the optimization-based methods. This difference is particularly noticeable in the cases without perturbations, meaning that the NRM can be the preferable choice in terms of accuracy when the input data does not present noise or artifacts. One of the advantages of using optimization methods is the treatment of the tracking constraint equations as objective functions to be minimized during the dynamic analysis. This approach ensures that only the topological constraints are treated as equality constraints, thus preventing violations of these constraints. All the sequential optimizations resulted in maximum constraint violations for the topological constraints on the order of 10 –10 and 10 –11 , within the number of iterations and function evaluation defined for the analysis and considering a constraint tolerance of 10 –10 . It should be noted that for some positions of the five-bar linkage model, the CNO with the interior-point algorithm presented convergence problems due to the system matrix becoming close to singular or badly scaled, reaching the maximum number of iterations. Despite reducing the average time of the analysis, the use of a parallel approach to solve the optimization problem resulted in a noisy pattern. This outcome can be explained by using the initial position as the optimization initial guess for all time frames. The method converged for similar positions that comply with the optimization equality constraints, but do not represent the global minimal point. Hence, the kinematic positions should be post-processed to avoid unrealistic accelerations, as these could impact the dynamic outcomes. An alternative solution to mitigate this issue is to first perform an analysis using the NRM and WLS and then use the results from this analysis as initial guesses for the optimization problem. Another approach is to subdivide the entire problem into as many sub-problems as the available threads, solving the kinematic analysis sequentially within each sub-problem. However, while this strategy addresses the issue of noisy positions within each sub-problem, it may introduce the challenge of ensuring data continuity at the boundaries between sub-problems. Hence, considering the advantages and limitations of the tested methods, optimization algorithms should be the preferred choice when analysis accuracy is critical, while the NRM with WLS should be employed when computation effort is an important factor, or the input data does not present noise or perturbations. It is worth mentioning that a 3-second analysis period sampled at 100 Hz was used in this study, as this duration exceeds the typical cycle times observed in human motion. For example, a full gait cycle at normal cadence lasts approximately 1.1 s [49–51]. While this duration might be considered short for simulation-based analyses [52], the present study focuses exclusively on the inverse dynamic analysis of (bio)mechanical systems, which does not require integration steps. Therefore, as long as the initial guess for each time frame is close to the actual solution, both the numerical and optimization-based methods are expected to converge without numerical problems. In the approach adopted here, the output from one time frame serves as the initial guess for the subsequent frame, further ensuring proximity to the solution and promoting fast convergence. This assumption may not hold for very fast motions; however, in such cases, an increased acquisition frequency is typically required to accurately capture the dynamics of the movement [53], thereby mitigating this issue. A relevant feature of the MC formulation with FCC-GRB is that it largely preserves the principles of the classical FCC-GRB approach. Therefore, the advantages reported for the classical formulation are still maintained. However, the inclusion of angular coordinates in the vector of the generalized coordinates implies that the system is no longer exclusively defined using Cartesian coordinates. Consequently, kinematic constraint equations incorporating nonlinear terms dependent on angular variables are generated. Since these variables do not directly define the rigid bodies, they do not impact the definition of the rigid body topological constraints or the computation of other points or vectors of the model, as the definition of the transformation matrices C and V remains equal. Therefore, most of the kinematic constraints maintain their linear or quadratic nature, meaning that their contributions to the Jacobian matrix and right-hand side vectors of velocities and accelerations are still null, linear or quadratic with respect to the generalized coordinates. This independence of the rigid bodies’ definition from the angular generalized variables implies that the definition of the system mass matrices is the same as in the classical formulation, considering the elimination method described in section 2.6. Although forward dynamics is not the primary focus of this hybrid formulation or of the analyses presented in this study, the MC formulation with the FCC-GRB framework has the potential to be adapted for the forward dynamics analysis of mechanical systems. While the increased dimensionality of the system is expected to raise computational demands, no additional numerical challenges are anticipated beyond those commonly encountered in traditional FCC-GRB simulations, such as singular configurations, redundant constraints, ill-conditioned Jacobian and mass matrices, or integration stability problems [29,41,54]. Future work also aims to explore extending the MC with FCC-GRB to incorporate flexible bodies or adopting alternative solving strategies based on different optimization algorithms or artificial intelligence approaches [37,38,55–58]. The similarities between the classical FCC-GRB formulation and the FCC-GRB with MC in terms of implementation and modeling procedure imply that the straightforwardness associated with the classical approach remains applicable in the second case. The potential applications identified for the classical formulation, such as teaching multibody dynamic topics in higher education, are still valid [28,29]. Although the MC variation includes angular variables, it does not require a high level of expertise in 3D rotation parametrization. This characteristic aligns with the advantages often highlighted when utilizing the natural coordinates formulation in advanced educational contexts [35]. The intrinsic properties of the FCC-GRB formulation also make it particularly suitable for the biomechanics of motion field. Since S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 22 the points and vectors that compose the generic rigid body have a direct relation with the inertial information found in anthropometrical databases, the modeling process becomes straightforward. The incorporation of the MC simplifies the inverse dynamic analysis of biomechanical models, as the reference points and vectors can be defined in the measure that they directly represent the markers used in the experimental acquisition of movement. Moreover, the rigid body vectors can be easily aligned with the anatomical joint axes, directly defining the rotation directions. If the generalized angular drivers are defined to match the convention established for the joint angles ([3,4]), the methodology can directly provide the joint angles. Finally, as highlighted in the introduction, experimental data acquired using traditional marker-based optoelectronic motion capture systems are prone to STA, which, due to their frequency characteristics and unpredictable nature, are difficult to automatically detect and attenuate using conventional smoothing techniques. By fitting the biomechanical multibody model, typically defined using measurements taken during static acquisitions with the markers placed at the correct anatomical landmarks, the MC helps mitigate STA effects, thereby improving both kinematic and dynamic analyses. This advantage is particularly significant, as both dynamic and muscle analysis outcomes are highly sensitive to perturbations in kinematic data [15,59,60]. Hence, in addition to simplifying the preand post-processing steps, the MC method enhances the accuracy of kinematic data, ultimately improving the reliability of biomechanical analysis. It is also important to mention that the use of the MC approach increases the complexity of the system, particularly when solved from a forward dynamics perspective, due to the increase in the number of coordinates and kinematic constraint equations. When applied in an inverse kinematic context, the increase of the system complexity has less impact on the computational performance of the analysis, making it more suitable for this type of study. In fact, the kinematic analysis with NRM and WLS for the five-bar linkage took approximately 30 ms per frame to converge, indicating its potential for real-time applications. This time could be improved by using a reduced definition of the model or by using a faster compiled programming language instead of an interpreted one. In certain scenarios, employing an MC formulation can offer additional modeling advantages when applied from a forward dynamics or predictive standpoint [35]. The explicit use of angular variables can simplify the control of specific systems that rely on angular inputs and facilitate the modeling of the action of external torque actuators, such as passive and active exoskeletons and wheelchairs [61–65]. 6. Conclusions This work explores the use of an MC formulation with FCC-GRB in the inverse kinematic and dynamic analysis of spatial multibody systems. The major steps required to fully implement the formulation, both with fully-defined and reduced rigid bodies, are described in detail. The formulation is validated in the analysis of two models with different levels of complexity. The accuracy and computational differences between using a numerical method based on the NRM and optimization algorithms are analyzed to evaluate the applicability of the formulation within different contexts. The results indicate that in terms of accuracy, optimization-based methods ensure the kinematic consistency of the system, albeit the increase in computational effort. The use of an NRM with WLS allows for the reduction of the computational cost, as it enables to control the violation of constraints at a topological level by defining appropriate weights. From an implementation point of view, the use of an MC formulation with FCC-GRB offers the main advantage of simultaneously computing the kinematically consistent position and drivers associated with the joint angular DoFs. This procedure is achieved by treating the joint displacement angles as generalized coordinates of the system and by minimizing the distance between a set of points and vectors of the model and their respective reference data. As a result, there is no need for preor post-processing steps to obtain all kinematic outcomes, avoiding also the errors associated with the use of angular drivers, such as the propagation of experimental errors along the kinematic chain. Despite requiring the explicit use of angular variables, the adoption of MC does not require complex parametrizations of 3D angle rotations, maintaining the simplicity of implementation typically attributed to the natural coordinates and classical FCC-GRB formulation. Therefore, it is the authors’ belief that this formulation remains suitable not only for teaching multibody-related topics but also for other areas that require kinematic analyses of complex systems with noisy experimental data, such as biomechanics of movement. CRediT authorship contribution statement S´ ergio B. Gonçalves: Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization. Ivo Roupa: Writing – review & editing, Validation, Conceptualization. Paulo Flores: Writing – review & editing, Supervision. Miguel Tavares da Silva: Writing – review & editing, Supervision, Project administration, Conceptualization. Declaration of competing interest The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: S´ ergio B. Goncalves, Ivo Roupa, Paulo Flores, and Miguel Tavares da Silva reports financial support was provided by Fundaç˜ ao para a Ciencia e a Tecnologia (FCT). If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 23 Acknowledgments The authors acknowledge Fundaç˜ ao para a Ciˆ encia e a Tecnologia (FCT) for its financial support via the projects LAETA Base Funding (DOI: 10.54499/UIDB/50022/2020), LAETA Programmatic Funding (DOI: 10.54499/UIDP/50022/2020), UIDB/04436/ 2020 and UIDP/04436/2020, and Portuguese Recovery and Resilience Program (PRR) for its financial support via IAPMEI/ANI/FCT under Agenda C645022399-00000057 (eGamesLab). Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.mechmachtheory.2025. 106080. Data availability No data was used for the research described in the article. References [1] M.P.T. Silva, J.A.C. Ambr´ osio, Kinematic data consistency in the inverse dynamic analysis of biomechanical systems, Multibody Syst. Dyn. 8 (2002) 219–239, https://doi.org/10.1023/A:1019545530737. [2] S. Ausejo, A. Suescun, J. Celigüeta, X. Wang, Robust Human motion reconstruction in the presence of missing markers and the absence of markers for some body segments, SAE Technical Papers (2006). https://doi.org/10.4271/2006-01-2321. [3] G. Wu, S. Siegler, P. Allard, C. Kirtley, A. Leardini, D. Rosenbaum, M. Whittle, D.D. D’Lima, L. Cristofolini, H. Witte, ISB recommendation on definitions of joint coordinate system of various joints for the reporting of human joint motion—Part I: ankle, hip, and spine, J. Biomech. 35 (2002) 543–548. [4] G. Wu, F.C.T. der Helm, H.E.J.D. Veeger, M. Makhsous, P. Van Roy, C. Anglin, J. Nagels, A.R. Karduna, K. McQuade, X. Wang, ISB recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion—Part II: shoulder, elbow, wrist and hand, J. Biomech. 38 (2005) 981–992. [5] S.B. Gonçalves, S.B.C. Lama, M.T. da Silva, Three decades of gait index development: a comparative review of clinical and research gait indices, Clin. Biomech. 96 (2022) 105682, https://doi.org/10.1016/J.CLINBIOMECH.2022.105682. [6] R. Stagni, S. Fantozzi, A. Cappello, A. Leardini, Quantification of soft tissue artefact in motion analysis by combining 3D fluoroscopy and stereophotogrammetry: a study on two subjects, Clin. Biomech. (Bristol,. Avon) 20 (2005) 320–329, https://doi.org/10.1016/j.clinbiomech.2004.11.012. [7] M. Begon, M.S. Andersen, R. Dumas, Multibody kinematics optimization for the estimation of upper and lower limb Human joint kinematics: a systematized methodological review, J. Biomech. Eng. 140 (2018) 030801, https://doi.org/10.1115/1.4038741. [8] A. Peters, B. Galna, M. Sangeux, M. Morris, R. Baker, Quantification of soft tissue artifact in lower limb human motion analysis: a systematic review, Gait Posture 31 (2010) 1–8, https://doi.org/10.1016/j.gaitpost.2009.09.004. [9] A. Leardini, A. Chiari, U. Della Croce, A. Cappozzo, Human movement analysis using stereophotogrammetry part 3. Soft tissue artifact assessment and compensation, Gait Posture 21 (2005) 212–225, https://doi.org/10.1016/j.gaitpost.2004.05.002. [10] J. Cl´ ement, R. Dumas, N. Hagemeister, J.A. de Guise, Soft tissue artifact compensation in knee kinematics by multi-body optimization: performance of subjectspecific knee joint models, J. Biomech. 48 (2015) 3796–3802, https://doi.org/10.1016/J.JBIOMECH.2015.09.040. [11] V. Radhakrishnan, M. Robinson, N.M. Fiorentino, S.B. Patil, A. Pelah, Reducing soft tissue artefacts through projection of markers and microwave imaging: an exploratory study, Sci. Rep. 15 (2025) 1–23, https://doi.org/10.1038/S41598-025-89586-W. [12] S. Guan, H.A. Gray, F. Keynejad, M.G. Pandy, Mobile biplane X-ray imaging system for measuring 3D dynamic joint motion during overground gait, IEEE Trans. Med. ImAging 35 (2016) 326–336, https://doi.org/10.1109/TMI.2015.2473168. [13] M. Schulze, F. Trautwein, T. Vordemvenne, M. Raschke, F. Heuer, A method to perform spinal motion analysis from functional X-ray images, J. Biomech. 44 (2011) 1740–1746, https://doi.org/10.1016/J.JBIOMECH.2011.03.040. [14] R. Dumas, V. Camomilla, T. Bonci, L. Cheze, A. Cappozzo, A qualitative analysis of soft tissue artefact during running, Comput. Methods Biomech. Biomed. Eng. 17 (2014) 124–125, https://doi.org/10.1080/10255842.2014.931518. [15] M.T. Silva, J.A.C. Ambr´ osio, Sensitivity of the results produced by the inverse dynamic analysis of a human stride to perturbed input data, Gait Posture 19 (2004) 35–49, https://doi.org/10.1016/S0966-6362(03)00013-4. [16] W. Schiehlen, Research trends in multibody system dynamics, Multibody Syst. Dyn. 18 (2007) 3–13, https://doi.org/10.1007/S11044-007-9064-4. [17] J.G. de Jal´ on, E. Bayo, Kinematic and Dynamic Simulation of Multibody Systems: the Real-Time Challenge, Springer Verlag, New York, 1994. [18] S.L. Delp, F.C. Anderson, A.S. Arnold, P. Loan, A. Habib, C.T. John, E. Guendelman, D.G. Thelen, OpenSim: open-source software to create and analyze dynamic simulations of movement, IEEE Trans. Biomed. Eng. 54 (2007) 1940–1950. [19] R. Dumas, L. Cheze, Soft tissue artifact compensation by linear 3D interpolation and approximation methods, J. Biomech. 42 (2009) 2214–2217, https://doi. org/10.1016/J.JBIOMECH.2009.06.006. [20] T. Ryu, H.S. Choi, M.K. Chung, Soft tissue artifact compensation using displacement dependency between anatomical landmarks and skin markers – a preliminary study, Int. J. Ind. Ergon. 39 (2009) 152–158, https://doi.org/10.1016/J.ERGON.2008.05.005. [21] A. Cappello, R. Stagni, S. Fantozzi, A. Leardini, Soft tissue artifact compensation in knee kinematics by double anatomical landmark calibration: performance of a novel method during selected motor tasks, IEEE Trans. Biomed. Eng. 52 (2005) 992–998, https://doi.org/10.1109/TBME.2005.846728. [22] M.S. Andersen, M. Damsgaard, B. MacWilliams, J. Rasmussen, A computationally efficient optimisation-based method for parameter identification of kinematically determinate and over-determinate biomechanical systems, Comput. Methods Biomech. Biomed. Eng. 13 (2010) 171–183, https://doi.org/ 10.1080/10255840903067080. [23] M.S. Andersen, M. Damsgaard, J. Rasmussen, Kinematic analysis of over-determinate biomechanical systems, Comput. Methods Biomech. Biomed. Eng. 12 (2009) 371–384, https://doi.org/10.1080/10255840802459412. [24] P. Nikravesh, Computer-aided Analysis of Mechanical Systems, 1st edition, 07632, Prentice Hall, New Jersey, 1988. [25] E.J. Haug, Computer aided Kinematics and Dynamics of Mechanical Systems, 1, Basic Methods, Allyn & Bacon, Inc., 1989. [26] A. Shabana, Dynamics of Multibody Systems, Cambridge University Press, 2020. [27] J.G. De Jal´ on, J. Unda, A. Avello, Natural coordinates for the computer analysis of multibody systems, Comput. Methods Appl. Mech. Eng. 56 (1986) 309–327. [28] I. Roupa, S.B. Gonçalves, M.T. da Silva, Kinematics and dynamics of planar multibody systems with fully Cartesian coordinates and a generic rigid body, Mech. Mach. Theory 180 (2023) 105–134, https://doi.org/10.1016/j.mechmachtheory.2022.105134. S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 24 [29] S.B. Gonçalves, I. Roupa, P. Flores, M. Tavares da Silva, Kinematic and dynamic analysis of spatial multibody systems with fully Cartesian coordinates and a generic rigid body, Mech. Mach. Theory 209 (2025) 1–35, https://doi.org/10.1016/j.mechmachtheory.2025.105955. [30] G. Gim, P.E. Nikravesh, Joint coordinate method for analysis and design of multibody systems: part 1. System equations, KSME J. 7 (1993) 14–25, https://doi. org/10.1007/BF02953141. [31] A. Seth, M. Sherman, P. Eastman, S. Delp, Minimal formulation of joint motion for biomechanisms, Nonlinear Dyn. 62 (2010) 291–303, https://doi.org/ 10.1007/s11071-010-9717-3. [32] M.T. Gameiro, P. Silva, Modelaç˜ ao e Simulaç˜ ao Sistem´ atica Em Coordenadas Cartesianas Totais de Sistemas Multicorpo, Actas Do Congresso de M´ etodos Num´ ericos Em Engenharia, Porto, Portugal, 2007, p. 2007. Junho 13-15. [33] J.A. Ambr´ osio, M. Tavares da Silva, A biomechanical multibody model with a detailed locomotion muscle apparatus. Advances in Computational Multibody Systems, Springer, Dordrecht, Netherlands, 2005, pp. 155–184. [34] M. Saura, J. Cuadrado, D. Dopico, A.I. Celdran, Computational kinematics of multibody systems: the advantages of a topological method based on its kinematic structure, in: Proceedings of ECCOMAS Multibody Dynamics, 2013, p. 1. [35] J.G. Jal´ on, Twenty-five years of natural coordinates, Multibody Syst. Dyn. 18 (2007) 15–33. [36] I. Roupa, M.R. da Silva, F. Marques, S.B. Gonçalves, P. Flores, M.T. da Silva, On the modeling of biomechanical systems for Human movement analysis: a narrative review, Arch. Comput. Methods Eng. 29 (2022) 4915–4958, https://doi.org/10.1007/s11831-022-09757-0. [37] N. Song, H. Peng, X. Guo, Sym-ML: a symplectic machine learning framework for stable dynamic prediction of mechanical system, Mech. Mach. Theory 206 (2025) 105934, https://doi.org/10.1016/J.MECHMACHTHEORY.2025.105934. [38] A. Hashemi, G. Orzechowski, A. Mikkola, J. McPhee, Multibody dynamics and control using machine learning, Multibody Syst. Dyn. 58 (2023) 397–431, https://doi.org/10.1007/S11044-023-09884-X. [39] S. Akram, Q.U. Ann, Newton raphson method, Int. J. Sci. Eng. Res. 6 (2015) 1748–1752. [40] T.J. Ypma, Historical development of the Newton-Raphson method, SIAM Rev. 37 (1995) 531–551, https://doi.org/10.1137/1037125. [41] J. García de Jal´ on, M.D. Guti´ errez-L´ opez, Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraint forces, Multibody Syst. Dyn. 30 (2013) 311–341, https://doi.org/10.1007/s11044-013-9358-7. [42] Y. Wang, Gauss–Newton method, Wiley Interdiscip. Rev. Comput. Stat. 4 (2012) 415–420, https://doi.org/10.1002/WICS.1202. [43] F.W. Gembicki, Vector Optimization for Control with Performance and Parameter Sensitivity Indices, Ph.D. Thesis, Case Western Reserve Univ, 1974. [44] A. Ancillao, E. Aertbeli¨ en, J. De Schutter, Effect of the soft tissue artifact on marker measurements and on the calculation of the helical axis of the knee during a gait cycle: a study on the CAMS-Knee data set, Hum. Mov. Sci. 80 (2021) 102866, https://doi.org/10.1016/J.HUMOV.2021.102866. [45] N.M. Fiorentino, P.R. Atkins, M.J. Kutschke, J.M. Goebel, K.B. Foreman, A.E. Anderson, Soft tissue artifact causes significant errors in the calculation of joint angles and range of motion at the hip, Gait Posture 55 (2017) 184, https://doi.org/10.1016/J.GAITPOST.2017.03.033. [46] D.A. Winter, H.G. Sidwall, D.A. Hobson, Measurement and reduction of noise in kinematics of locomotion, J. Biomech. 7 (1974) 157–159, https://doi.org/ 10.1016/0021-9290(74)90056-6. [47] F. Marques, I. Roupa, M.T. Silva, P. Flores, H.M. Lankarani, Examination and comparison of different methods to model closed loop kinematic chains using lagrangian formulation with cut joint, clearance joint constraint and elastic joint approaches, Mech. Mach. Theory 160 (2021) 104294, https://doi.org/ 10.1016/j.mechmachtheory.2021.104294. [48] R. Soram, S. Roy, S.R. Singh, M. Khomdram, S. Yaikhom, S. Takhellambam, On the rate of convergence of Newton-Raphson method, Int. J. Eng. Sci. (IJES) 2 (2013) 5–12. [49] J. Perry, J.R. Davids, Gait analysis: normal and pathological function, J. Pediatr. Orthopaed. 12 (1992) 815. [50] D.A. Winter, The Biomechanics and Motor Control of Human Gait: Normal, Elderly and Pathological, 2nd Edition, University of Waterloo Press, Waterloo, Ontario, Canada, 1991. [51] S.B. Gonçalves, M.T. Silva, J.M. Martins, M.C. Neves, Advanced Computer Methods for Pathological and Non-Pathological Human Movement Analysis, EUROMECH Colloquium 511 on Biomechanics of Human Motion, 2011. [52] M. Gonz´ alez, F. Gonz´ alez, A. Luaces, J. Cuadrado, A collaborative benchmarking framework for multibody system dynamics, Eng. Comput. 26 (2010) 1–9, https://doi.org/10.1007/s00366-009-0139-0. [53] F. Fallahtafti, S.R. Wurdeman, J.M. Yentes, Sampling rate influences the regularity analysis of temporal domain measures of walking more than spatial domain measures, Gait Posture 88 (2021) 216–220, https://doi.org/10.1016/J.GAITPOST.2021.05.031. [54] P. Flores, M. Machado, E. Seabra, M. Tavares da Silva, A parametric study on the Baumgarte Stabilization method for forward dynamics of constrained multibody systems, J. Comput. Nonlinear Dyn. 6 (2011) 011019, https://doi.org/10.1115/1.4002338. [55] J.A.C. Ambr´ osio, M.A. Neto, R.P. Leal, Optimization of a complex flexible multibody systems with composite materials, Multibody Syst. Dyn. 18 (2007) 117–144, https://doi.org/10.1007/S11044-007-9086-Y. [56] V. Gufler, E. Wehrle, A. Zw¨ olfer, A review of flexible multibody dynamics for gradient-based design optimization, Multibody Syst. Dyn. 53 (4) (2021) 379–409, https://doi.org/10.1007/S11044-021-09802-Z. [57] Y. He, J. McPhee, A design methodology for mechatronic vehicles: application of multidisciplinary optimization, multibody dynamics and genetic algorithms, Veh. Syst. Dyn. 43 (2005) 697–733, https://doi.org/10.1080/00423110500151077. [58] N. Song, M. Wang, X. Wang, H. Peng, A novel machine learning method for real-time dynamic analysis of tensegrity flexible multibody systems, Nonlinear Dyn. (2025) 1–28, https://doi.org/10.1007/S11071-025-11152-W. [59] S.B. Gonçalves, M.R. da Silva, F. Marques, P. Flores, M.T. da Silva, Validation of skeletal muscle models in multibody dynamics: a collaborative collection of benchmark cases, Submitted to. Multibody System Dynamics, 2025. [60] C. Redl, M. Gfoehler, M.G. Pandy, Sensitivity of muscle force estimates to variations in muscle–tendon properties, Hum. Mov. Sci. 26 (2007) 306–319, https:// doi.org/10.1016/J.HUMOV.2007.01.008. [61] W.J. Jaimes, J.F. Mantilla, S.A. Salinas, H.J. Navarro, Modeling and Simulation of a Lower Limb Exoskeleton with Computed Torque Control for Gait Rehabilitation, Pan American Health Care Exchanges, PAHCE, 2021, https://doi.org/10.1109/GMEPE/PAHCE50215.2021.9434854, 2021-May. [62] K.A. Inkol, J. McPhee, Assessing control of fixed-support balance recovery in wearable lower-limb exoskeletons using multibody dynamic modelling, in: Proceedings of the IEEE RAS and EMBS International Conference on Biomedical Robotics and Biomechatronics, 2020, pp. 54–60, https://doi.org/10.1109/ BIOROB49111.2020.9224430, 2020-November. [63] S. Doung, U. Wasiwitono, Multibody dynamics modeling and control of wheelchair balancing system, in: Proceedings - 2021 International Seminar on Intelligent Technology and Its Application: intelligent Systems for the New Normal Era 2021, ISITIA, 2021, pp. 123–128, https://doi.org/10.1109/ ISITIA52817.2021.9502215. [64] L.P. Quinto, S.B. Gonçalves, M.T. Silva, Design of a passive exoskeleton to support sit-to-stand movement: a 2D model for the dynamic analysis of motion. Biosystems and Biorobotics, Springer, Cham, 2019, pp. 299–303, https://doi.org/10.1007/978-3-030-01887-0_57. [65] M. Harant, M.B. N¨ af, K. Mombaur, Multibody dynamics and optimal control for optimizing spinal exoskeleton design and support, Multibody Syst. Dyn. 57 (2023) 389–411, https://doi.org/10.1007/S11044-023-09877-W. S.B. Gonçalves et al. Mechanism and Machine Theory 214 (2025) 106080 25