A general framework for modeling and dynamic simulation of multibody systems using factor graphs

In this paper, we present a novel general framework grounded in the factor graph theory to solve kinematic and dynamic problems for multi-body systems. Although the motion of multi-body systems is considered to be a well-studied problem and various methods have been proposed for its solution, a unified approach providing an intuitive interpretation is still pursued. We describe how to build factor graphs to model and simulate multibody systems using both, independent and dependent coordinates. Then, batch optimization or a fixed-lag-smoother can be applied to solve the underlying optimization problem that results in a highly-sparse nonlinear minimization problem. The proposed framework has been tested in extensive simulations and validated against a commercial multibody software. We release a reference implementation as an open-source C++ library, based on the GTSAM framework, a well-known estimation library. Simulations of forward and inverse dynamics are presented, showing comparable accuracy with classical approaches. The proposed factor graph-based framework has the potential to be integrated into applications related with motion estimation and parameter identification of complex mechanical systems, ranging from mechanisms to vehicles, or robot manipulators.


Introduction
Dynamic simulation of multibody systems mainly refer to inverse and forward problems. Inverse dynamics deals with the determination of the driving forces that generate a given motion, as well as the constraint reaction forces. The solution to the inverse dynamics problem can be obtained for example using the Newton-Euler (N-E) method [1], that results in efficient recursive algorithms, e.g., Rigid Body Dynamics Library (RBDL) [2] and similar methods [3]. Conversely, forward dynamics involves the motion estimation of a multi-body system over time under the given applied forces and initial conditions. Therefore, in a direct dynamic problem, the motion is the result of the force system that generates it. From a mathematical perspective, forward dynamics is computationally intensive as it entails the integration of a system of nonlinear ordinary differential equations. The most common formulations that deal with forward dynamic are the Composite-Rigid-Body Algorithm (CRBA) [4] and the Articulated-Body Algorithm (ABA) [5]. However, the above mentioned methods are not suited for systems with closed-loops like parallel robots, for which more complex and expensive procedures are required to solve both inverse and forward dynamics, as described for example in [6]. Therefore, a single algorithm that can solve all types of dynamics problems has not been fully established. One notable effort is the work by Rodriguez [7], who proposed a unified approach based on Kalman filtering and on the concept of smoothing to solve dynamics problems. This approach has been further extended to solve the forward dynamic problem of closed kinematic chains [8]. Another effort in defining a unifying framework has been given by [9], who analyzed various algorithms for serial chain dynamics. In [10], an attempt is proposed to unify the CRBA and ABA derivation as two elimination methods to solve forward dynamics. To the best of our knowledge, factor graphs have not been applied to solve the general multibody kinematics and dynamics in the existing literature. The closest works are the preprints [11,12], where factor graphs are indeed employed to solve kinematic and dynamic problems, although their proposed graph structure is applicable to openloop robot arms only.
The present work advocates factor graphs as a unifying graphical language in which to express the classical kinematics and dynamics algorithms, with the additional potential to develop novel and advanced state estimators. The main contributions of this work are: • a factor graph-based representation of dynamics problems, which is a insightful visualization of the underlying sparse constraints between all involved variables, • a unified method, which can solve inverse and forward dynamics for either open or closed kinematic chains, • a flexible framework that can be expressed and solved for both dependent and independent coordinates.
Also, note that despite the implementation described in this manuscript is focused on planar mechanisms, it is perfectly suitable to spatial systems without changes at the level of factor graphs. Additionally, our approach allows more powerful and flexible schemes for state and parameter estimation to be implemented in contrast with standard methods based on Kalman filtering [13,14,15,16,17]. However, such applications are left for future extensions of this work to keep the present manuscript focused on the key ideas on how to apply factor graphs to multibody motion problems.
The proposed approach draws on the formalism of graphical models, a powerful tool borrowing concepts from statistics and graph theory [18,19]. By addressing the multibody simulation problem from the perspective of the variable structure, graphical models allow us creating efficient estimators for any combination of observed and hidden variables, effectively unifying the problems of kinematic and forward dynamic analysis (predicting or estimating the trajectory of a MB system), inverse dynamics, and parameter identification (e.g., inertial properties of the bodies involved, external disturbance forces, friction in the joints, etc.). All those problems end up to be formulated as a sparse nonlinear cost function built from a library of reusable "building blocks" (the factors) on which efficient solvers can then be applied. The framework is suitable for either offline batch analysis or online real-time operation that represents another clear advantage of the proposed approach.
The rest of the paper is structured as follows. Sections 2 and 3 first provide the required background on multibody dynamics and graphical models, respectively. Next, section 4 presents a methodology for the application of Bayesian networks to multibody dynamics problems, whereas section 5 particularizes such networks as factor graphs for a number of practical problems. Individual factors used in those graphs are described in detail in section 6. Numerical examples are provided in section 7, and we end sketching some conclusions in section 8.

Review of Multibody Dynamics
In this section, fundamentals of multibody system motion analysis are briefly recalled, whereas the interested reader can refer for more details to the wide related Literature, e.g. [20]. A multibody system is an assembly of two or more bodies (or elements) constrained to each other to fulfill a given motion law. In many practical applications, these elements may be considered rigid and, throughout this paper, we will work under this assumption even though the proposed framework may be further extended to include body flexibility. One of the key decisions to take when modeling a MBS is selecting the set of generalized coordinates that will be used to represent it. Using independent coordinates z allows one to deal with the lowest number of parameters, i.e. the number of DOFs of the system. A second choice is to adopt dependent coordinates q, in a number larger than that of DOFs but able to describe all multibody system points univocally. When dependent coordinates are used, the corresponding set of constraint equations must be included as well for a complete system analysis. There are different kinds of dependent coordinates [20]: Relative Coordinates, References Point Coordinates, Natural Coordinates, and a combination of the previous ones (Mixed Coordinates). Natural coordinates, mixed with relative coordinates where needed, will be assumed in this work.
In the remainder of this Section, the MBS motion equations are developed and expressed in terms of both dependent and independent coordinates.

Dependent coordinates formulation
For any given MBS with f dofs, the use of n dependent coordinates expressed by the vector q requires m = n − f constraint equations, which form the following set of equations Thus, it is assumed that there are at least as many equations as there are unknown coordinates. To solve the kinematic problem, the time derivative of Eq. (1) is required, one time for velocity analysis and two times for acceleration analysis, leading to the following set of equations whereq is the vector of dependent velocities, Φ q ∈ R m×n the Jacobian matrix of Eq. (1) respect to q, and Φ t the time derivative of constraint equations that is equal to the null vector for time-independent constraints. To solve the dynamic problem, external forces and inertia forces need to be considered. From the classical Newton's law, one obtains the following set of differential equations that expresses the force equilibrium equations where λ is the vector of Lagrange multipliers, M is the system mass matrix and vector Q contains the generalized external forces. Since Eq. (4) is a system of n equations in n + m variables, by adding the m equations (3) one obtains the following system Normally [20,21], this equation is solved for the extended vector of unknowns that includes both,q and λ. In our framework, we are specifically interested in the generalized accelerationsq. Therefore, by applying the block matrix inversion lemma (see 8) to Eq. (5) and keeping the first row only, it leads to the following equation of motion expressed in terms of dependent coordinates: where we defined the auxiliary term Γ (q) as the m × m square symmetric matrix introduced for convenience in subsequent derivations.

Independent coordinates formulation
Independent coordinates z, ensure the minimum number of variables, i.e. the number of DOFs. However, multibody systems can be more difficult to analyze with respect to dependent ones. First, the matrix R is introduced aṡ where the columns of R are f linearly independent vectors that constitute a basis of the nullspace of Φ q . Then, we express the independent velocitiesż as the projection of the dependent velocitiesq on the rows of a constant matrix B that we assume satisfying the condition of having f rows linearly independent from each other and from the m rows of Φ qż = Bq By combining Eq. (2) with Eq. (9), one obtains The columns of matrix R are the partial velocities with respect to the generalized coordinates z, and the term Sb represents the partial velocities with respect to time. By differentiating Eq. (9) and by taking into account Eq. (3), one gets From Eq. (4), pre-multiplied by R R Mq = R Q (13) and by substitutingq in Eq. (12) Introducing the reduced mass matrix M = R MR and force vector Q = R (Q − MSc), the following equation of motion in terms of independent coordinates can be finally obtained

Background on graphical models
A factor graph is a particular type of probabilistic graphical model that can be used to describe the structure of an estimation problem [22]. Factor graphs have found applications in many fields, for example in robot perception [23], information theory [24], signal processing [25], and in other areas of robotics, mostly SLAM [26] and computer vision [27], state estimation in legged robots [28], and kinematic motion planning [29].

Dynamic Bayesian Networks
It is instructive to start our discussion by considering a Dynamic Bayesian Network (DBN), a type of graphical model where variables are represented as nodes and directed edges stand for causal relationships [30,22]. The existence of directed edges in a DBN allows us to encode an expert's knowledge in causality relationships between all involved variables in the graph.
The rules to convert a DBN into the kind of graphs we are actually using in this work, factor graphs (FGs), are discussed in section 5. A relevant point regarding such conversion is that one single DBN can be mapped into different FGs depending on which subset of the original DBN variables are known data and which unknowns are required to be found. Let us remark that whereas a DBN displays all variables as nodes, in a FG only unknown variables are variable nodes.

Factor graphs
Factor graphs (FGs) are bipartite graphs comprising two kinds of nodes: variable nodes, and factor nodes. Variables are the unknown data to be estimated, and factors represent any constraint (a cost function to be minimized) or relationship between variables. A crucial aspect of a FG is its sparsity: each factor node is only connected to the variables that appear in its cost function. Sparse graph optimization algorithms exist with computational cost almost linear with the number of edges in a FG, as opposed to, for example, the cubic cost of a naive implementation of Kalman Filtering [31]. Incremental (e.g. [32]) and sliding-window (e.g. [33]) approaches exist as well, enabling the efficient estimation of problems with thousands to millions of variables [34,35,36].
Once a FG is formulated for a given problem, estimating the most-likely value of all the unknowns becomes a numerical nonlinear minimization problem, for which very efficient algorithms exist. As a probabilistic estimator, these optimal values can be assigned an uncertainty. In general, retrieving the covariance Σ of the estimated variables involves inverting the Hessian of the linearized problem evaluated at the optimal solution, such that Σ = (J ΛJ) −1 with J the sparse Jacobian of all constraints (factors), and Λ the weight matrix representing our level of certainty about each constraint. The matrix inversion operation is, in general, very costly since it is cubic with the problem dimension. However, optimized methods exist for the kind of sparse problem at hand, see for example, [37, §B.4] or [32].

DBN for multibody dynamics
Since DBNs allow a human expert to specify the causality relationships between variables, we propose the two DBNs in Figure 1 as the underlying structure of a generic multi-body system that is the basis for our proposed framework. Note the use of discrete time steps t, exactly as simulations are commonly done following numerical integration schemes. A set of observations or measurements o are available from sensors installed on the system. Note, however, Model for independent coordinates formulation that the graphical model formalism is flexible enough to allow each sensor to be available at a different or even sampling rate. When using this graphical model to achieve state estimation, sensors may provide partial information on the system dynamics, which will be then fused with the rest of graph constraints over the system trajectory to reach the most-likely values of the estimated variables. In non state-estimation problems, the set of observations o may be empty. Additionally, system parameters such as masses, stiffness or friction coefficients (constant or variable over time) could be added as additional nodes, although this possibility is left out of the scope of this manuscript for the sake of conciseness. From Figure 1, some observations can be drawn: • A central role in the DBN for independent coordinates (Figure 1(b)) is played by the independent coordinates z,ż,z as they represent the degrees of freedom (d.o.f) that govern the evolution of the MBS. This role is assumed by the dependent coordinate variables q,q,q in Figure 1(a).
• The branch variable C, required uniquely in the independent-coordinate formulation, allows us to disambiguate between, e.g. the two possible configurations for a four-bar linkage if we are only given the information about the crank angle. This variable could be part of the unknowns to be estimated, as already done in [38], although that work did not use the more powerful FG representation proposed here.
• Typically, nodes in a DBN are depicted as shaded or unshaded depending of whether they are "hidden" or "observed" variables, respectively, e.g. see [19,38]. The former are estimated from the latter. Since this work focuses on FGs instead, where observed variables are subsumed into factor nodes and unknowns become variable nodes [23], different FGs will be generated from the same DBNs in Figure 1 for different multibody simulation problems, hence it is not convenient to establish such a distinction at the DBN level yet.
• Observations o (sensor readings) are a function of (all, or a subset of) dependent coordinates. Typical observations that may be useful in a MBD problem are measurements obtained from gyroscopes, accelerometers, and load cells.
• External forces Q t act by means of modifying accelerationsq t orz t , according to the system dynamics, expressed by Eq. (6) or Eq. (15), respectively.
Once the model of the dynamics MBS has been expressed as a graphical model, the time evolution of the system can be obtained by converting the graph into a maximum-a-posteriori (MAP) estimation problem. To this aim, one could write down the joint probability distribution p(φ) of all the involved variables φ = {q 0:t ,q 0:t ,q 0:t , o 0:t , Q 0:t } for time steps 0 to t exploiting the conditional probability encoded by the DBN (refer to e.g. [19,22]), which for dependent coordinates (i.e. Figure 1(a)) becomes: where each conditional probability term in Eq. (16), taking negative logarithms and up to an irrelevant proportionality constant, becomes a factor in this alternative form of an overall cost function c(φ) to be minimized: A fundamental feature of graphical models is how they enable factorizing functions of all problem variables such as c(φ) ∝ − log p(φ) into the product of a large number of smaller functions (in terms of the number of variables involved in each expression) called factors in Eq. (17), which are discussed individually in section 6. This is what keeps the estimation problem tractable even for hundreds of thousands of variables, something intractable for estimators in the family of the Kalman filter not exploiting the sparsity of the problem structure. Note that the goal of an estimator is searching for the optimal set of unknowns φ that maximize the likelihood of all observed data according to the model, that is: Depending on the division of the DBN variables into observed and unknowns, we would arrive at different factorizations since factors are considered to be functions of the unknowns only. For example, if all variables in one given cost factor in Eq. (17) are known data it just evaluates to a constant which can be absorbed by the proportionality relationship c(φ) ∝ − log p(φ) and hence can be ignored. A graphical representation of the remaining relevant factors leads us to factor graphs themselves.

Factor graphs for multibody dynamics
The conversion from DBN to FG is known to be achievable as follows, without the need to explicitly writing down probabilistic equations like Eq. (16)-Eq. (17): Every Bayes net [DBN] node splits in both a variable node and a factor node in the corresponding factor graph. The factor is connected to the variable node, as well as the variable nodes corresponding to the parent nodes in the Bayes net. If some nodes in the Bayes net are evidence nodes, i.e., they are given as known variables, we omit the corresponding variable nodes: the known variable simply becomes a fixed parameter in the corresponding factor. [23, p. 12] Applying these rules to the DBN reflecting the structure of variables involved in the dynamic simulation of multibody systems with dependent coordinates in Figure 1(a) we derive next the FGs for a couple of practical problems depending on which are the known and unknown variables. Factors will be only briefly discussed here regarding their purpose and meaning, whereas their detailed implementations are described in section 6.

FG for forward dynamics in dependent coordinates
In forward dynamics simulation, we are given a known initial state of a mechanical system (position q 0 and velocity q 0 ), their geometric and inertial properties, and the external forces that act on it (Q i , i = 0, ..., N). For simplicity, we assume no sensors are installed in the system since this problem instance does not need them, but predicting sensor outputs could be also possible adding the corresponding nodes and factors.
The resulting FG when using dependent coordinates is shown in Figure 2 and it involves the following factor nodes (or plain factors): Figure 2: Factor graphs for the forward dynamic simulation problem using dependent generalized coordinates. Circle nodes are problem unknowns, solid square nodes are factors. See discussion in section 5.1.
• f prior : Factors imposing a given a priori knowledge about the attached variables, e.g. initial conditions. Prior factors can be defined for both, position and velocity.
• f d pc : Factor for position constraints in dependent coordinates. It ensures the fulfillment of mechanical holonomic constraints by keeping the state vector q on the manifold Φ(q, t) = 0, hence this factor is repeated for each position node q i . It could be omitted for q 0 if the enforced pose imposed by the prior factor is known to be a correct mechanism position complying with Φ(·) = 0; alternatively, the prior factor can be made to affect a minimum number of variables in q 0 (the number of degrees of freedom), leaving f d pc in charge of recalculating the rest of the generalized coordinates. This factor, on its own, solves the so-called position problem [20] in multibody dynamics.
• f d vc : Factor for velocity constraints in dependent coordinates, enforcing generalized velocitiesq to remain on the manifold Φ qq + Φ t = 0.
• f d dyn : Factor for the dynamics equation of motion: it links external forces (known data in this FG, hence not reflected as variable nodes) with the instantaneous accelerationq. It also depends on q andq since acceleration is always a function of them too. Note in the graph how acceleration for each timestep depends on position and velocity of the same timestep only.
• f ti : Trapezoidal numerical integrator factor is used twice per timestep: to integrate velocities into positions, and accelerations into velocities. Implicit integrators as the Trapezoidal one are often employed in MBS simulations for their stability. However, the Euler integrator version is also devised (see section 6), since explicit integrators are commonly used in real-time applications as well.

FG for forward dynamics in independent coordinates
Alternatively, one could devise a FG implementation for the forward dynamics problem when using independent coordinates, leading to the graph depicted in Figure 3. In this case, the following factors are used: • f prior : In this case, they are used to impose both, an initial known dynamic state (z 0 ,ż 0 ) and approximated initial values for the full vector of dependent coordinates q 0 . The latter may be required to solve ambiguities in closed kinematic topologies, e.g. a four bar linkage, where knowing the minimum set of independent variables still leaves more than one feasible kinematic configuration.
• f eq : An equality factor, used to impose a soft constraint between state vectors consecutive in time. The rationale behind this factor is to impose a soft constraint, which can be easily violated (its weight is small in comparison to all other factors) but still provides a solid starting point for nonlinear numerical solvers to exploit the fact that mechanism coordinates cannot vary abruptly between consecutive time steps. This factor is particularly important to solve ambiguities in mechanisms with more than one branch, exactly as argued by the end of the former paragraph. • f i dyn : The independent-coordinate version of f d dyn discussed above, it implements the dynamics equation of motion.
• f i pc , f i vc : Like their dependent-coordinate counterparts in the former section, these factors keep the position and velocity vector within their corresponding constraint manifolds. Note how, in this case, the factors not only affect q andq coordinates, but their independent-coordinate versions z andż too, respectively. These factors actually correspond to the so-called position problem and velocity problem in multibody mechanics [20].

FG for inverse dynamics in dependent coordinates
Another problem that can be formalized as a FG is inverse dynamics, as shown in Figure 4. The problem consists of specifying a desired trajectory over time for the set of degrees of freedom of a mechanism, then solving for the required forces and torques to generate such trajectory. The following factors are used: • f prior : Here, different prior factors are used to define initial known values for the dynamic state (q 0 ,q 0 ) and to enforce a value of zero in all components of the generalized force vectors Q i where it is known that no external force is acting. In other words, the latter factors are required to leave only part of Q i as an unknown, e.g. for the degree of freedom where a motor is exerting a torque.
• f d idyn : Similar to the dynamics equation factor f d dyn discussed above, but with the force Q i as an additional variable. Note that f d dyn (see section 5.2) also used a value for Q i but it was treated as a known constant, whereas for f d idyn the Q i are unknowns and as such the factor provides a Jacobian of the error term with respect to them as well.
From all discussed problems so far, inverse dynamics is the hardest to solve for nonlinear iterative solvers from initial values that are far from the optimum. Therefore, it requires a proper initialization strategy to enable numerical nonlinear solvers to cope with it effectively, as discussed in the experiments section later on.

Factors definition
In a factor graph, factors establish the relations between variables. This section provides an insightful practical guide to those factors MBS problems are made of. Each factor defines an error e(x) between predicted and measured data. To apply nonlinear optimization algorithms (e.g. Gauss-Newton or Levenberg-Marquardt), the Jacobian of such error with respect to all involved variables x is required. DBN variables whose values are known do not become FG nodes and hence are constant parameters of factors. The goal of the optimization is to search for the best values of all variables x taking the weighted sum of error functions (one per factor) as close to zero as possible: where e i (x) 2 Σ is a form of whitening already integrating the probabilistic noise model (or weight) of each factor. The most common model, used in the present work, is the assumption of additive zero mean Gaussian noise n ∼ N(0, Σ) with covariance matrix Σ. Taking the negative logarithm of the corresponding probability density can be shown to give us the nonlinear least-squares problem in Eq. (19), where with information matrix Λ = Σ −1 . Larger information values in Λ (or smaller variances in Σ) imply that the factor must be considered with a higher weight during the optimization in comparison to other factors with smaller information values (or larger variances). Note that each component of multidimensional x vectors may have a different weight for a given single factor, e.g. as proposed in the priors over external forces in section 5.3. The interested reader is deferred to [23] for a more in-depth discussion on this topic. In the following, note the use of the superscripts f d · and f i · for factors with differentiated implementations for dependent and independent-coordinates, respectively.

f prior : Prior factor
Prior factors f prior (x) over a variable x are the only unary factor used in our proposed graphs. They represent a priori belief (hence the name) about the state of a given variable.
Since the problems addressed in this manuscript use variables that are either (a) state vectors with generalized coordinates of planar mechanisms, or (b) their velocities or (c) their accelerations, all variables can be treated as vectors in the group of real numbers, hence the error function e(·) of this factor for a vector of dimensionality d has a simple form: f prior error function: Variables: x (21b) Fixed parameters: x 0 (21c) Jacobian: where x 0 is the "desired" value for the variable x.

f ei : Euler integrator factor
The graphs proposed in former sections work over temporal sequences of variables, sampled at a fixed sample rate f s = 1/∆t. Imposing the continuity of ordinary differential equations in discrete time can be done by means of numerical integration, of which the Euler integrator is the simplest instance.
Given two consecutive samples for time steps t and t + 1 of a given variable x and its time derivativeẋ at time t, the Euler integrator factor f ei (shown in Figure 5(a)) can be defined as: Variables: Fixed parameters: ∆t (22c) Jacobians: where the state space is assumed to be R d .

f ti : Trapezoidal integrator factor
Another well-known numerical integration scheme follows the Trapezoidal Rule. Given two consecutive samples for time steps t and t + 1, both for a given variable x and its time derivativeẋ, the Trapezoidal integrator factor f ti , shown in Figure 5(b), can be defined as: Variables: Fixed parameters: ∆t (23c) Jacobians: As mentioned in section 5.1, this integrator is preferred for its better accuracy in comparison to the Euler method.

f d pc : Factor for position constraints in dependent coordinates
This factor, depicted in Figure 5(c), ensures the fulfillment of mechanical holonomic constraints of Eq. (1). As explained in Section 2.1, modeling a mechanism in dependent coordinates leads to a number of constraint equations largest or equal to the number of generalized coordinates.
This factor has the following general form: Fixed parameters: (constant distances, angles, etc.) (24c) Jacobians: where both Φ and Φ q can be automatically built out of pre-designed blocks (see Appendix III), according to a formal description of the topology of the mechanical system and the specific joints connecting each pair of adjacent bodies. In particular, since each physical constraint becomes one or more entries in the Φ(q) vector, each such constraint fully determines the corresponding rows in the Jacobian Φ q (q), easing its automated construction from the elemental expressions exposed in Appendix III.
6.5. f d vc : Factor for velocity constraints in dependent coordinates Velocity constraint factors further improve the quality of MBS simulation. These kind of factors are modeled by Eq. (2) and depend on both, q andq, as illustrated Figure 5(d). Its definition is as follows: Fixed parameters: (constant distances, angles, etc.) (25c) Jacobians: where we assumed that no constraint depends explicitly on time, hence Φ t (q t ) can be safely neglected. Again, each row in the Jacobians comes from exactly one constraint in the definition of the mechanical system, with most common cases described in Appendix III.
6.6. f d dyn , f d idyn and f i dyn , f i idyn : Factors for dynamics These factors minimize the error between the actual acceleration estimates (q andz) and the corresponding equations of motion Eq. (6) and Eq. (15) for dependent and independent coordinates, respectively. To deal both forward and inverse dynamics, this factor connects the generalized positions, velocities, and accelerations with the forces (Q t ) for each single time step t. In the case of forward dynamics, both the variable Q t and the edge connecting it to the factor are dropped, becoming Q t an internal known parameter to the factor. For this reason, two factors for forward and inverse dynamics must be defined independently:

• Forward dynamics
For dependent coordinates: Fixed parameters: (masses, external forces) (26c) Jacobians: and for independent coordinates: Fixed parameters: (masses, external forces) (27c) Jacobians: • Inverse dynamics For dependent coordinates: Fixed parameters: (masses) (28c) Jacobians: and for independent coordinates: Fixed parameters: (masses) (29c) Jacobians: Refer to Eq. (6) and Eq. (15) for the evaluation of the error functions above. Numerical Jacobians have been used in our implementation for Eq. (28d) and Eq. (29d), since exact closed-form Jacobians are so complex than the computational advantage of having closed-form expressions. In this particular case is not clear and deserves future analysis. Figure 3 shows three sets of factors linking dependent coordinates (position, velocities, and accelerations) with their corresponding degrees of freedom, that is, with their counterparts in independent coordinates. The first one of these factor is f i pc , in charge of imposing the simultaneous fulfillment of: (a) position constraints in Φ(q), and (b) that independent coordinates z are integrated into its corresponding positions within q. Therefore, its error function and Jacobians read: Jacobians:

f i pc : Factor for position constraints in independent coordinates
where {idxs} = {y 1 , y 2 , ..., y d } stands for the fixed sequence of indices of each z coordinate inside the n-vector q, and the coefficients I i, j of the binary matrix I idxs are defined as 1 if j = y i , or as 0 otherwise, where i = 1, ..., d and j = 1, ..., n.
6.8. f i vc : Factor for velocity constraints in independent coordinates This factor ensures the fulfillment of the velocity constraints for independent coordinates: Fixed parameters: {idxs}, constant distances, angles, etc. Jacobians: where the tensor-vector product Φ qqqt results in the 2-D matrixΦ q (see Eq. (40a)) which can be systematically built row by row from a description of the mechanism from the building blocks described in Appendix III.
6.9. f i ac : Factor for acceleration constraints in independent coordinates This factor ensures the fulfillment of the acceleration constraints for independent coordinates: Fixed parameters: (constant distances, angles, etc.) (32c) Jacobians: where the tensor-vector products (Φ q ) q (q t )q t , Φ qqqt , and 2Φ qqqt = 2Φ q can be systematically built as described in Appendix III.

Test case validation
A four-bar linkage is employed to exemplify the implementation and the performance of the proposed FG-based framework. A scheme of the mechanism is shown in Figure 6(a), where the two revolute joints, P 1 and P 2 , can move in the motion plane, whereas A and D are fixed. Geometric and inertial properties are summarized in Figure 6(b). The motion of the mechanism has been simulated using the commercial MBS environment Adams/View from MSC.Software with a fixed-step integrator at 1/10 of the time step used in the factor graph, hence its results can be considered the ground-truth for comparison purposes.
A companion open-source reference C++ implementation 1 of the proposed factor graph approach has been released along with the paper. Note that all factor closed-form Jacobians have been thoroughly validated by means of unit tests against numerical finite differences. θ P 1 : ( x 1 , y 1 ) Figure 6: (a) The planar mechanism used in the numeric simulations described in section 7. (b) Table of mechanism properties (length and mass) for each link. Inertias I i are given with respect to the center of mass of each link, placed at its center.

Forward dynamic simulation
Free motion of the four-bar mechanism under the sole effect of gravity (g = 9.8 m/s 2 ) has been simulated, with an initial crank angle θ (refer to Figure 6(a)) of zero and null initial velocities.
The FG in Figure 2 using dependent coordinates has been implemented, and then solved numerically in an incremental fashion using a fixed-lag smoother, when only the factors of the last N w time steps (the window length) up to time t are considered, with discrete time t ranging from 0 up to 5 seconds with steps of dt = 1 ms. A Levenberg-Marquardt nonlinear iterative solver algorithm has been selected for these runs, with a maximum number of iterations of 15. Note, however, that virtually all calls to the optimizer required only 3 to 5 iterations until convergence. Figure 7(a)-(b) show the time history of the dependent generalized coordinates of points P 1 and P 2 , and their corresponding velocities, respectively, obtained from both, the proposed FG-based approach and from MSC Adams. Both series accurately coincide, hence a more detailed comparison is in order in Figure 7(c)-(d), where the differences in q andq are reported between the proposed FG-based method and the ground truth. Similar good agreement (omitted for brevity's sake) is obtained when adopting the independent coordinates FG, as explained in Figure 3, thus verifying our two proposed schemes. The noise models used in all factors involved are summarized in Table 2. Covariance absolute values are not relevant, only their relative weights with respect to each other.
In order to quantify the accuracy of both schemes (dependent and independent-based FGs) and to explore the effect of the fixed-lag window length parameter N w , the root mean square error (RMSE) of the overall q trajectories are shown in Table 1 for window size varying from 2 to 10 time steps. Two main conclusions can be drawn from these data: (a) the dependent-coordinates scheme seems to provide a higher accuracy, and (b) using window lengths larger than a few (3-4) time steps do not further increase the accuracy of the estimation. This latter point can be explained by the fact that the present simulations focus on forward dynamics only, not estimation from noisy sensor measurements, a case where it should be expected that larger windows lead to better results. This topic, however, falls out of the scope of the present paper and it should be studied in the future.    Table 2: Noise models (covariance matrices) employed in the forward dynamics simulations. Most covariances represent isotropic models (those including the identity matrix I, or 0 for perfectly-known values), with only a few of them having different weights for different coordinates (those defined with a diag(· · · ), a diagonal matrix with the given diagonal values).

Inverse dynamic simulation
Next, the FG proposed in Figure 4 for inverse dynamic calculation is put at test with the same four-bar mechanism described above. Here, the crank angle θ trajectory over time is specified as an arbitrary smooth curve, and the goal is to solve for the unknown torque at the crank (point A) that would be required to balance gravity and inertial forces as accurately following the prescribed path.
In this case, we use a Levenberg-Marquardt nonlinear solver on the FG in Figure 4 in batch mode, as opposite to filtering or fixed-lag smoother mode. However, the position-problem is strongly nonlinear, preventing solvers to find the global optimal of the whole FG in Figure 4 for arbitrary initial values for all variables. Instead, we propose attacking the inverse dynamics FG in the following four stages (running the batch optimizer once after each stage), which ensure that the global optimum is always easy to find by the nonlinear solver: • Stage 1: Include q variables, prior factors for the desired trajectory ( f prior ), position-constraint factors ( f d pc ), and soft equality factors ( f eq ).
• Stage 3: Add velocity-constraint factors ( f d vc ). • Stage 4: Add generalized force variables Q and inverse dynamics factors ( f d idyn ). Numerical results for this approach are shown in Figure 8, where the trajectories prescribed and actually-followed by the mechanism crank are depicted in Figure 8(a) for both the proposed approach and the commercial software MSC ADAMS. Both methods successfully solve the inverse dynamics problem by finding the torque sequence in Figure 8(b). In this case, our method leads to an excellent stable following error (less than 1 milliradian), whereas MSC's error seems to slowly grow over time, as shown in Figure 8(c).

Conclusions and future works
This work has settled the theoretical bases for formulating kinematics and dynamics problems from computational mechanics in the form of factor graphs, and demonstrated the validity of our reference open-source implementation. The results showed the practical feasibility of the proposed approach and its accuracy in comparison with a commercial MBS software based on classical approaches. One additional advantage of our formalism lies at its excellent suitability to design state observers and parameter identification systems for arbitrary mechanical systems, ranging from robots to vehicles. This aspect will be part of future works.

Appendix I. Block matrix inversion lemma
This lemma implies that:

Appendix II. Tensor notation definitions and useful expressions
Order of a tensor: Number of indices to index all its components, i.e. Scalars, vectors, and matrices are 0th, 1st, and 2nd-order tensors, respectively.
Einstein convention: Superscript indices imply reading elements up-down, subscript indices imply reading leftright. Summation works over the repeated index or indices. Example with a common matrix product: Conventional notation: C ik = j A i j B jk Einstein convention: Transpose rule: We have that A i j = (A ) j i . Jacobian: ∂x ∂y , assuming x and y are a row and a column vector, respectively, becomes: The standard derivative chain rule: Applied to vector variables x and z, using an intermediary variable y, can be put in tensor notation as:   Chain rule for derivatives of matrices with respect to a vector: It gives us the third-order tensor: where it becomes clear that we must sum over the l indices, that is, converting the Einstein notation above into a conventional summation: Product of a third-order tensor and a vector: Given a third-order tensor T = T i j k and a vector v = v k , its product reads: Appendix III. Equations for planar multibody systems In the following, we summarize the equations implemented in our open-source library to model typical mechanical constraints of planar mechanisms. These equations are the building blocks of the errors and Jacobian matrices associated with the factors described in sections 6.4-6.
5. An open-source implementation of all these equations is provided in the C++ toolkit mentioned in section 7.
Since each constraint contributes only a few scalar components to Φ (typically only one for planar mechanisms), the following uses the superscript Φ k to reflect the k-th row in the corresponding vector or matrix, with k being arbitrarily determined by the order in which constraints are enumerated when defining a mechanism. Furthermore, for the sake of notation clarity, q should be understood in the following sets of equations as the subset of the actual q that contains only those generalized coordinates that are directly involved in the constraint, i.e. the provided Jacobians should be expanded with columns of zeros as needed at the position of non-relevant coordinates.
In addition to the specific formulas for each constraint, the following terms appear in factor Jacobians but are generic so they apply to all constraints:

Constant distance
This constraint is typically used to model rigid solids. Given the coordinates of any pair of points i and j of a rigid body, such that (x i , y i ) ∈ q and (x j , y j ) ∈ q, this constraint ensures a fixed distance L between them. Each such constraint contributes with the following rows in the different vectors and matrices defining the multibody system: Fixed pinned slider constraint This constraints model a point P(x, y) enforced to move along a line defined by two fixed points A(x A , y A ) and B(x B , y B ). By exploiting the properties of the similar triangles, one obtains: Mobile pinned slider constraint A more general version of the former constraint, where the point P(x, y) is a now moving along a line defined by two mobile points A(x i , y i ) and B(x j , y j ). In this case: Absolute orientation of planar link Another special kind of constraints is that one needed when defining relative coordinates, such as the absolute angle of a link, e.g. a mechanism input crank, or a vehicle wheel position. The relative angle discussed in this section is called absolute since it is defined as the angle θ between a body comprising a pair of points (x i , y i )-(x j , y j ), at a fixed distance L, and the horizontal axis of the global frame of reference.