*2.1. The Multibody Equations of Motion*

The mathematical description of the system dynamics can be derived by means of the Lagrange's equations for constrained mechanical systems [18]:

$$\frac{d}{dt}\left(\frac{\partial \mathcal{L}(\dot{q}, q\_{\prime}\lambda)}{\partial[\dot{q}, \lambda]}\right) - \frac{\partial \mathcal{L}(\dot{q}, q\_{\prime}\lambda)}{\partial[q\_{\prime}\lambda]} = u^{\varepsilon}(q\_{\prime}u)\_{\prime} \tag{1}$$

with the Lagrangian defined as:

$$
\mathcal{L} = \mathcal{T} - \mathcal{V} - \phi(q)^T \lambda. \tag{2}
$$

L represents the Lagrangian functional, T the kinetic energy, V the potential energy, *φ*(*q*)*Tλ* the constraint contribution with the Lagrange multipliers *λ*, and *u<sup>e</sup>* is the vector of the external actions. MB models describe the dynamics of several rigid and/or flexible interacting bodies linked together through the definition of kinematic joints which are mathematically represented by the constraint equations *φ*(*q*) while *J* = *∂φ*(*q*)/*∂q* represent the Jacobian of the constraint equations. *<sup>q</sup>* <sup>∈</sup> <sup>R</sup>*nq* is the generalized coordinates vector, *<sup>λ</sup>* <sup>∈</sup> <sup>R</sup>*n<sup>λ</sup>* is known as Lagrange multipliers and *<sup>u</sup>* <sup>∈</sup> <sup>R</sup>*nu* is the input vector. Through the definition of the assembled body coordinates and the motion parametrization Equation (1) can be written in a residual form as a fully implicit real-valued non-linear function:

$$
\lg(\vec{q}, \vec{q}, q, \lambda, \mathfrak{u}) = 0.\tag{3}
$$

#### *2.2. The Differential-Algebraic form of the EOMs*

A set of natural coordinates *qn* <sup>∈</sup> <sup>R</sup>*nqn* was proposed in [18], where redundant degrees of freedom are employed to define the system coordinates of the assembled bodies. Moreover, including the motion parameterization employed in the Flexible Natural Coordinate Formulation (FNCF) [19] allows deriving a constant singular mass matrix *<sup>M</sup>* <sup>∈</sup> <sup>R</sup>*nqn*×*nqn* . Assuming this formulation, Equation (3) can be written in the so-called index-3 form:

$$\begin{cases} \mathbf{g}\_1 = M\_n \ddot{\mathbf{q}}\_n + f\_{nl}(\dot{\mathbf{q}}\_{n\prime} q\_{n\prime} \boldsymbol{\mu}) + f^T \boldsymbol{\lambda} = \mathbf{0}\_q\\ \mathbf{g}\_2 = \boldsymbol{\phi}(q\_n) = \mathbf{0}\_{\boldsymbol{\lambda}} \end{cases} \quad . \tag{4}$$

Here, *fnl* <sup>∈</sup> <sup>R</sup>*nqn* is the non-linear generalized force vector expressed as:

$$f\_{\rm nl}(\phi\_{\rm nl}, q\_{\rm nl}, \boldsymbol{u}) = f\_{\rm \upsilon}(\phi\_{\rm n}, q\_{\rm n}) + f\_{\rm int}(\phi\_{\rm n}, q\_{\rm n}) + f\_{\rm ext}(\phi\_{\rm n}, q\_{\rm n}, \boldsymbol{u}).\tag{5}$$

*fv* represents the quadratic velocity vector related to the gyroscopic forces of the bodies, which is zero for FNCF formulation. *fint* is the internal force vectors which accounts for the elastic energy stored by deformable bodies and if rigid bodies are assumed *fint* vanishes; *fext* is the external force vector and can be spilt in the sum of two contributions, the interaction forces among bodies *fb* (i.e., contact and friction forces) and the input forces *fu*. They can be summarized as follows:

$$f\_{\text{ext}}(\dot{q}\_{\text{n}}, q\_{\text{n}}, \mathbf{u}) = f\_{\text{b}}(\dot{q}\_{\text{n}}, q\_{\text{n}}) + f\_{\text{u}}(q\_{\text{n}}, \mathbf{u}). \tag{6}$$

Here, *fu* can be written as *fu*(*qn*, *u*) = *Ut*(*qn*)*u*, where *Ut* is tangent input matrix defined as:

$$
\Omega I\_t = \frac{\partial f\_u}{\partial u}.\tag{7}
$$

Due to the structure of the EOM, Equation (4), for the FNCF formulation, derivatives can be more readily obtained than for may alternative flexible multibody formulations. Therefore, the above mentioned coordinates definition and motion parameterization will be considered in this work. For the sake of brevity, we omit the subscript *n* referred to the natural coordinate formulation for the remainder of this manuscript.

Despite the computational advantages of the above mentioned MB approach, the methodologies that will be introduced in the next sections can be easily extended to alternative MB formulation, such as the floating-frame of reference component mode synthesis approach or the generalized component mode synthesis [5,20].

The I-DAEs form of Equation (4) are generally not suitable for estimation algorithms such as the Kalman Filter family, since these have been designed to handle E-ODEs type of equations.

In the next section, we present a new methodology to directly linearize the I-DAEs starting from its discrete form but without employing any explicit constraint elimination technique.
