2.1.1. Double-Step Semi-Recursive Formulation

In the semi-recursive formulation, a body *i* is defined by the set of six Cartesian velocities as **Z***<sup>i</sup>* = **r**˙ *T <sup>i</sup> <sup>ω</sup><sup>T</sup> i T* and six Cartesian accelerations as **Z**˙ *<sup>i</sup>* = **r**¨*T <sup>i</sup> <sup>ω</sup>*˙ *<sup>T</sup> i T* for a complete description [61,62]. Here, **r**˙*i*, **r**¨*i*, *ω<sup>i</sup>* and *ω*˙ *<sup>i</sup>* are velocities, accelerations, angular velocities, and angular accelerations of the body, respectively. In the relative coordinate system, the position of *nb* bodies in a system can be described by using joint coordinates as **z** = *z*<sup>1</sup> *z*<sup>2</sup> .... *znb T* [61,62]. The absolute velocity **Z** and the absolute acceleration **Z**˙ of the system bodies can be mapped in terms of the relative joint velocity vector **z**˙ and the relative joint acceleration vector **z**¨ by using the velocity transformation matrix as follows [61,62]:

$$\begin{aligned} \mathbf{Z} &= \mathbf{T} \mathbf{R}\_d \dot{\mathbf{z}} \\ \dot{\mathbf{Z}} &= \mathbf{T} \mathbf{R}\_d \ddot{\mathbf{z}} + \mathbf{T} \dot{\mathbf{R}}\_d \dot{\mathbf{z}} \end{aligned} \tag{1}$$

where **<sup>T</sup>** <sup>∈</sup> <sup>R</sup>6*nb*×6*nb* is the path matrix that demonstrates the topology of the system, and **<sup>R</sup>***<sup>d</sup>* <sup>∈</sup> <sup>R</sup>6*nb*×*nb* is a block diagonal matrix. The path matrix **<sup>T</sup>** is a lower triangular matrix and contains entries of 6 × 6 (**I**6) unit matrices representing bodies between the body under observation and the root of the system [61]. In Equation (1), the block diagonal matrix **R***<sup>d</sup>* and the product **R**˙ *<sup>d</sup>***z**˙ can be computed with the joint-dependent element of the velocity transformation matrix **<sup>b</sup>***<sup>i</sup>* <sup>∈</sup> <sup>R</sup>6×<sup>1</sup> and the joint-dependent element of the acceleration transformation vector **<sup>d</sup>***<sup>i</sup>* <sup>∈</sup> <sup>R</sup>6×1, respectively [61,62]. The semi-recursive formulation is described hereafter, but the interested reader is referred to [61,62] for further details of **T**, **b***i*, and **<sup>d</sup>***i*. The composite mass matrix **<sup>M</sup>***<sup>i</sup>* <sup>∈</sup> <sup>R</sup>6×<sup>6</sup> and the composite force vector **<sup>Q</sup>***<sup>i</sup>* <sup>∈</sup> <sup>R</sup>6×<sup>1</sup> of the *i*th body can be described using absolute coordinates as [61,62]:

$$\overline{\mathbf{M}}\_{i} = \begin{bmatrix} m\_{i}\mathbf{I}\_{3} & -m\_{i}\overline{\mathbf{g}}\_{i} \\ m\_{i}\overline{\mathbf{g}}\_{i} & \mathbf{J}\_{i} - m\_{i}\overline{\mathbf{g}}\_{i}\overline{\mathbf{g}}\_{i} \end{bmatrix} \tag{2}$$

$$\overline{\mathbf{Q}}\_{i} = \begin{bmatrix} \overline{\mathbf{f}}\_{i} - \tilde{\omega}\_{i}(\tilde{\omega}\_{i} m\_{i} \mathbf{g}\_{i})\\ \mathbf{n}\_{i} - \tilde{\omega}\_{i} \mathbf{l}\_{i} \omega\_{i} + \tilde{\mathbf{g}}\_{i} (\mathbf{f}\_{i} - \tilde{\omega}\_{i} (\tilde{\omega}\_{i} m\_{i} \mathbf{g}\_{i})) \end{bmatrix}' \tag{3}$$

where *mi* is the mass of the *<sup>i</sup>* body, **<sup>f</sup>***<sup>i</sup>* <sup>∈</sup> <sup>R</sup>3×<sup>1</sup> is a vector of external forces, *<sup>ω</sup>*˜ *<sup>i</sup>* is the skew-symmetric matrix of the angular velocity vector, **<sup>n</sup>***<sup>i</sup>* <sup>∈</sup> <sup>R</sup>3×<sup>1</sup> is the vector of external moments, **<sup>I</sup>**<sup>3</sup> is a 3 <sup>×</sup> 3 unit matrix, and **<sup>g</sup>**˜*<sup>i</sup>* <sup>∈</sup> <sup>R</sup>3×<sup>3</sup> is the skew-symmetric matrix of the position vector of the centre of mass of the body in the global coordinate system. In Equation (2), **J***<sup>i</sup>* is the inertia tensor of the *i*th body, which can be computed as described in [61]. Applying the principle of virtual work and using Equation (1) yields the equation of motion for an open-loop system in the simplified form [61,62]:

$$\mathbf{R}\_d \prescript{T}{}{\mathbf{T}}^T \overline{\mathbf{M}} \mathbf{T} \mathbf{R}\_d \dot{\mathbf{z}} = \mathbf{R}\_d \prescript{T}{}{\mathbf{T}}^T (\overline{\mathbf{Q}} - \overline{\mathbf{M}} \mathbf{T} \dot{\mathbf{R}}\_d \dot{\mathbf{z}}),\tag{4}$$

where **<sup>M</sup>** <sup>∈</sup> <sup>R</sup>6*nb*×6*nb* is the block diagonal matrix consisting of the composite mass matrices of the bodies. The force vector **<sup>Q</sup>** <sup>∈</sup> <sup>R</sup>6*nb*×<sup>1</sup> is the column matrix of composite forces. To incorporate closed-loop systems, the double-step semi-recursive formulation is used in this study [62]. In this method, a set of *m* constraint equations **Φ**(**z**) = **0** are introduced for closure of an open-loop system. This method employs Gaussian elimination with a full pivoting approach to identify the independent and dependent columns of the Jacobian matrix **Φ<sup>z</sup>** [62–64]. Through this formulation, relative joint-independent coordinates can be used to define the dynamics of a system, i.e., the relative joint-dependent coordinates can be computed in terms of the relative joint-independent coordinates. Hence, this provides an appropriate option for the state and parameter estimation applications. The relative joint velocity vector ˙**z** can be described using the coordinate partitioning method [59]:

$$
\begin{bmatrix}
\dot{\mathbf{z}}^{\rm d} \\
\dot{\mathbf{z}}^{\rm i}
\end{bmatrix} = \begin{bmatrix}
\mathbf{I}
\end{bmatrix} \dot{\mathbf{z}}^{\rm i} \equiv \mathbf{R}\_{\mathbf{z}} \dot{\mathbf{z}}^{\rm i},\tag{5}
$$

where **<sup>z</sup>**˙ <sup>d</sup> <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* are the relative joint-dependent velocities, **<sup>z</sup>**˙ <sup>i</sup> <sup>∈</sup> <sup>R</sup>*nf* are the relative jointindependent velocities, **<sup>R</sup>**<sup>z</sup> <sup>∈</sup> <sup>R</sup>*nb*×*nf* is the velocity transformation matrix, **<sup>Φ</sup>**<sup>d</sup> **<sup>z</sup>** <sup>∈</sup> <sup>R</sup>*m*×*m*, and **Φ**<sup>i</sup> **<sup>z</sup>** <sup>∈</sup> <sup>R</sup>*m*×*nf* are the Jacobian matrices of the constraint equations with respect to the dependent and independent relative joint positions, respectively. In Equation (5), it is assumed that neither singular configurations nor redundant constraints exist; as a consequence, the inverse of matrix **Φ**<sup>d</sup> **<sup>z</sup>** can be obtained [55,59]. The relative joint acceleration vector can be expressed by differentiating Equation (5) [59]:

$$
\ddot{\mathbf{z}} = \mathbf{R}\_{\mathbf{z}} \ddot{\mathbf{z}}^{\mathrm{i}} + \mathbf{R}\_{\mathbf{z}} \dot{\mathbf{z}}^{\mathrm{i}},\tag{6}
$$

where **z**¨<sup>i</sup> are the relative joint-independent accelerations, and **R**˙ <sup>z</sup> is the derivative of the velocity transformation matrix. Substituting Equation (6) into Equation (4) results in an equation of motion for a closed-loop system in a simplified form [55,56,58,63,64]:

$$\mathbf{R}\_{\mathbf{z}}^T \mathbf{R}\_d^T \mathbf{T}^T \overline{\mathbf{M}} \mathbf{T} \mathbf{R}\_d \mathbf{R}\_{\mathbf{z}} \ddot{\mathbf{z}}^i = \mathbf{R}\_{\mathbf{z}}^T \mathbf{R}\_d^T \left(\mathbf{T}^T \overline{\mathbf{Q}} - \mathbf{T}^T \overline{\mathbf{M}} \mathbf{D}\right), \tag{7}$$

where **D** = **TR***<sup>d</sup>* ⎡ ⎣− **Φ**<sup>d</sup> **z** −1 **Φ**˙ **zz**˙ **0** ⎤ <sup>⎦</sup> <sup>+</sup> **TR**˙ *<sup>d</sup>***z**˙ represent the absolute accelerations when the vector **z**¨ is zero in Equation (6). Equation (7) can be further simplified using the accumulated mass matrix **M**<sup>Σ</sup> = **R***<sup>T</sup>* <sup>z</sup> **R***<sup>T</sup> <sup>d</sup>* **<sup>T</sup>***T***MTR***d***R**<sup>z</sup> and the accumulated force matrix **Q**<sup>Σ</sup> = **R***<sup>T</sup>* <sup>z</sup> **R***<sup>T</sup> d* **<sup>T</sup>***T***<sup>Q</sup>** <sup>−</sup> **<sup>T</sup>***T***MD** .

#### 2.1.2. Hydraulic Lumped Fluid Theory

The lumped fluid theory can be used to compute pressures within a hydraulic circuit [65]. Using this approach, a hydraulic circuit is assumed to be composed of discrete volumes. The pressures inside the volumes are equally distributed, with the acoustic waves having a negligible effect on the pressures [41,65]. In any hydraulic volume *Vh*, the differential pressure *p*˙ *<sup>h</sup>* can be computed [41,65] as

$$
\dot{p}\_h = \frac{k\_p + p\_h k\_0}{V\_h} \sum\_{h=1}^{n\_f} \mathbb{Q}\_{h\prime} \tag{8}
$$

where *Qh* is the sum of incoming and outgoing volume flow rates, *k*<sup>0</sup> is the flow gain, *kp* is the pressure flow coefficient, and *nf* is the total amount of volume flows. By employing a semi-empirical method, the volume flow rate *QR* through a throttle valve can be described as [66]

$$Q\_R = \mathbb{C}\_R \text{sgn}(\Delta p) \sqrt{|\,\Delta p\,\,|},\tag{9}$$

where Δ*p* is the pressure difference and *CR* = *CdAR* "2 *<sup>ρ</sup>* is the semi-empirical flow rate coefficient of the throttle valve. Here, *Cd* is the flow discharge coefficient and *ρ* is the fluid density. In Equation (9), *AR* is the cross-sectional area of the pressure relief valve. Similarly, the volume flow rate *Qd* through a directional control valve can be computed as [67,68]

$$Q\_d = \mathbb{C}\_v \iota \mathcal{U} \operatorname{sgn}(\Delta p) \sqrt{|\,\Delta p\,\,|},\tag{10}$$

where *Cv* is the semi-empiric flow rate coefficient, and *U* is the relative position of the spool. Equation (10) is complemented by the following first-order differential equation:

$$
\dot{\mathcal{U}} = \frac{\mathcal{U}\_{ref} - \mathcal{U}}{\tau},
\tag{11}
$$

where *Uref* is the reference voltage signal, and *τ* is the time constant. The incoming flow rate *Qin* and outgoing flow rate *Qout* in the hydraulic cylinder can be described as

$$Q\_{in} = \dot{s}A\_1, \qquad Q\_{out} = \dot{s}A\_2. \tag{12}$$

where *s*˙ is the piston velocity, and *A*<sup>1</sup> and *A*<sup>2</sup> are the areas on the piston and piston-rod side of the cylinder, respectively. The force *Fh* produced by the cylinder can be written as

$$F\_h = p\_1 A\_1 - p\_2 A\_2 - F\_{\mu\_{\prime}} \tag{13}$$

where *p*<sup>1</sup> and *p*<sup>2</sup> are, respectively, the pressure on the piston and piston-rod side, which can be calculated by using Equation (8). *Fμ* is the total friction force in the hydraulic cylinder caused by the hydraulic sealing. As proposed in [69], this friction force can be calculated by employing the Brown and McPhee model [70], which is valid for both positive and negative tangential velocity. The actuator velocity dependent friction force can be written in the vector form as

$$\mathbf{F}\_{\mu} = \left( F\_c \tanh\left(4 \frac{\|\mathbf{s}\|}{\mathbf{v}\_s}\right) + (F\_s - F\_c) \frac{\frac{\|\mathbf{s}\|}{\mathbf{v}\_s}}{\left(\frac{1}{4} \left(\frac{\|\mathbf{s}\|}{\mathbf{v}\_s}\right)^2 + \frac{3}{4}\right)^2} \right) \text{sgn}(\dot{\mathbf{s}}) + \sigma\_2 \dot{\mathbf{s}} \tanh(4), \qquad (14)$$

where *Fc* is the Coulomb friction, v*<sup>s</sup>* is the Stribeck velocity, *Fs* is the static friction, *σ*<sup>2</sup> is the coefficient of viscous friction, and **s**˙ is the actuator velocity vector.

#### 2.1.3. Monolithic Approach: Coupling MBS and Hydraulic Dynamic Systems

The MBS formulation can be combined with the fluid power system solver to form a unified set of non-linear differential equations in a monolithic approach:

$$\begin{aligned} \mathbf{M}^{\Sigma}(\mathbf{z})\ddot{\mathbf{z}}^{i} &= \mathbf{Q}^{\Sigma}(\mathbf{z}, \dot{\mathbf{z}}, \mathbf{p}) \\ \dot{\mathbf{p}} &= \mathbf{v}(\mathbf{z}, \dot{\mathbf{z}}, \mathbf{p}, \mathbf{y}, \mathsf{U}) \end{aligned} \qquad\qquad\qquad\qquad\qquad\tag{15}$$

where **p** is the pressure vector, and **y** is the vector of hydraulic parameters. Equation (15) is a set of non-linear equations that can be represented as **f**(**x**, *U*) = **0**. Here, the vector **x** = **z***<sup>T</sup>* **z**˙ *<sup>T</sup>* **p***<sup>T</sup>* **y***<sup>T</sup> T* . The solution of the non-linear equations described in Equation (15) is stiff. A stiff equation can be solved by using single-step implicit trapezoidal integration scheme [55,71–73]. In this integration scheme, the relative joint-independent positions and the pressures are initially predicted as **z**<sup>i</sup> *<sup>k</sup>*+<sup>1</sup> = **<sup>z</sup>**<sup>i</sup> *<sup>k</sup>* + **<sup>z</sup>**˙ <sup>i</sup> *<sup>k</sup>*Δ*<sup>k</sup>* <sup>+</sup> <sup>1</sup> 2 **z**¨i *<sup>k</sup>*Δ*k*<sup>2</sup> and **<sup>p</sup>***k*+<sup>1</sup> = **<sup>p</sup>***<sup>k</sup>* + **<sup>p</sup>**˙ *<sup>k</sup>*Δ*k*, respectively [73]. The derivatives of **z**<sup>i</sup> *<sup>k</sup>*+<sup>1</sup> and **p***k*+<sup>1</sup> can be predicted as

$$\begin{aligned} \dot{\mathbf{z}}\_{k+1}^{\dot{\mathbf{i}}} &= \frac{2}{\Delta k} \mathbf{z}\_{k+1}^{\dot{\mathbf{i}}} + \dot{\mathbf{z}}\_{k}^{\dot{\mathbf{i}}}\\ \dot{\mathbf{z}}\_{k+1}^{\dot{\mathbf{i}}} &= \frac{4}{\Delta k^2} \mathbf{z}\_{k+1}^{\dot{\mathbf{i}}} + \dot{\mathbf{z}}\_{k}^{\dot{\mathbf{i}}}\\ \dot{\mathbf{p}}\_{k+1} &= \frac{2}{\Delta k} \mathbf{p}\_{k+1} + \overset{\scriptstyle \vee}{\mathbf{p}}\_{k} \end{aligned} \tag{16}$$

where ˇ **z**˙ i *<sup>k</sup>* <sup>=</sup> <sup>−</sup>( <sup>2</sup> <sup>Δ</sup>*<sup>k</sup>* **<sup>z</sup>**<sup>i</sup> *<sup>k</sup>* + **<sup>z</sup>**˙ <sup>i</sup> *<sup>k</sup>*), <sup>ˇ</sup> **z**¨i *<sup>k</sup>* <sup>=</sup> <sup>−</sup>( <sup>4</sup> <sup>Δ</sup>*k*<sup>2</sup> **<sup>z</sup>**<sup>i</sup> *<sup>k</sup>* <sup>+</sup> <sup>4</sup> <sup>Δ</sup>*<sup>k</sup>* **<sup>z</sup>**˙ <sup>i</sup> *<sup>k</sup>* + **<sup>z</sup>**¨<sup>i</sup> *<sup>k</sup>*) and <sup>ˇ</sup> **<sup>p</sup>**˙ *<sup>k</sup>* <sup>=</sup> <sup>−</sup>( <sup>2</sup> <sup>Δ</sup>*k***p***<sup>k</sup>* + **p**˙ *<sup>k</sup>*). Note that the relative joint-dependent positions **z**<sup>d</sup> *<sup>k</sup>*+<sup>1</sup> are obtained from **<sup>z</sup>**<sup>i</sup> *<sup>k</sup>*+<sup>1</sup> and the previous step **<sup>z</sup>**<sup>d</sup> *k* by solving the position problem **Φ**(**z**) = 0 [59–61]. The non-linear constraint equations are solved iteratively with the Newton–Raphson method [59–61]. The derivatives of the relative joint-dependent positions **z**<sup>d</sup> *<sup>k</sup>*+<sup>1</sup> are computed from Equations (5) and (6) at the velocity and acceleration levels, respectively [55]. Substituting Equation (16) into Equation (15) leads to a set of dynamic equilibrium equations as **F**(*χk*+1) = **0** at the time step *k* + 1 as

$$\begin{aligned} \mathbf{M}^{\Sigma} \mathbf{z}\_{k+1}^{\mathrm{i}} - \frac{\Delta k^2}{4} \mathbf{Q}\_{k+1}^{\Sigma} + \frac{\Delta k^2}{4} \mathbf{M}^{\Sigma} \dot{\mathbf{z}}\_k^{\mathrm{i}} &= \mathbf{0} \\ \frac{\Delta k}{2} \mathbf{p}\_{k+1} - \frac{\Delta k^2}{4} \mathbf{v}\_{k+1} + \frac{\Delta k^2}{4} \dot{\mathbf{p}}\_k &= \mathbf{0} \end{aligned} \tag{17}$$

where *χk*+<sup>1</sup> = (**z**i )*T <sup>k</sup>*+<sup>1</sup> **<sup>p</sup>***<sup>T</sup> k*+1 *T* is unknown. The Newton–Raphson method is employed on the non-linear Equation (17) to iteratively compute the unknown variables [73,74]:

$$
\begin{bmatrix}
\frac{\partial \mathbb{F}(\boldsymbol{\chi})}{\partial \boldsymbol{\chi}}
\end{bmatrix}\_{k+1}^{(h)} \Delta \boldsymbol{\chi}\_{k+1}^{(h)} = -\left[\overline{\mathbf{F}}(\boldsymbol{\chi})\right]\_{k+1'}^{(h)}\tag{18}
$$

where 0 0 0 Δ**z**<sup>i</sup> *k*+1 0 0 <sup>0</sup> <sup>&</sup>lt; <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>10</sup> rad and<sup>0</sup> 0Δ**p***k*+<sup>1</sup> 0 <sup>0</sup> <sup>&</sup>lt; <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>2</sup> Pa are the error tolerances in the relative joint independent positions and pressures provided in the Newton–Raphson method. In Equation (18), **F**(*χ*) (*h*) *k*+1 is the residual vector, which can be computed as

$$\mathbb{E}\left[\mathbf{F}(\boldsymbol{\chi})\right]\_{k+1}^{(h)} = \frac{\Delta k^2}{4} \left[\mathbf{M}^{\Sigma}\mathbf{\bar{z}}^{\mathbf{i}} - \mathbf{Q}^{\Sigma}\right]\_{k+1}^{(h)}.\tag{19}$$

In Equation (18), *∂***F**(*χ*) *∂χ* (*h*) *<sup>k</sup>*+<sup>1</sup> is the tangent matrix, which can be numerically approximated at a point *χ*<sup>0</sup> by using a forward finite differences, as demonstrated in the literature [72,75]. Now, by computing Δ*χ*(*h*+1) *<sup>k</sup>*+<sup>1</sup> from Equation (18), the next iteration *<sup>χ</sup>*(*h*+1) *k*+1 can be calculated.

#### *2.2. Estimation Algorithm: ADEKF with a Curve-Fitting Method*

In this section, the ADEKF parameter estimation algorithm is introduced in the framework of a B-spline curve-fitting method. It is important to note, however, that the introduced procedure can be easily modified for applications of other curve-fitting methods, as mentioned in [48,49]. Parameter estimation through the ADEKF comprises prediction and update stages. At the prediction stage, in the case of the coupled multibody and hydraulic systems, the augmented state vector can be described as **x** = (**z**i )*<sup>T</sup>* (**z**˙ <sup>i</sup> )*<sup>T</sup>* **p***<sup>T</sup>* **y***<sup>T</sup> T* . At this step, **x**ˆ− *<sup>k</sup>* is calculated in the time integration of a dynamic model described as [46]

$$\hat{\mathbf{x}}\_{k}^{\prime -} = \mathbf{f}(\hat{\mathbf{x}}^{\prime +} \_{k-1}, \mathcal{U}\_{k})\_{\prime} \tag{20}$$

To account for unknown parameters, the proposed parameter estimation algorithm employs the curve-fitting method. Through this method, a B-spline curve is constructed with the knot vector **u** for non-uniform open splines [48,49] at the current time step:

$$\mathbf{C}(\mathbf{u}) = \sum\_{i=0}^{n} B^{i,d}(\mathbf{u}) \mathbf{N}^{i},\tag{21}$$

where n is the number of control points, *d* is the degree, *Bi*,*d*(**u**) are the *d*th order of B-spline basis functions, and **N***<sup>i</sup>* is the control point vector. The control point vector can be expressed in terms of the system parameters **y**. For instance, in the case of the characteristic curve, the control point vector can be written in terms of the spool position and semi-empiric flow rate coefficient as **N** = - *U*min *U*<sup>1</sup> ... *U*<sup>n</sup> *U*max *Cv*min *Cv*<sup>1</sup> ... *Cv*<sup>n</sup> *Cv*max . Here, *U*min, *U*1, and *U*<sup>n</sup> represent spool positions, and *Cv*min , *Cv*<sup>1</sup> , and *Cv*max are the semi-empiric flow rate coefficients of a hydraulic valve. *Bi*,*d*(**u**) can be defined by using the Cox–de Boor recursion formula [48,49]:

$$B^{i,0}(\mathbf{u}) = \begin{cases} 1 & \mathbf{u}^i \le \mathbf{u} < \mathbf{u}^{i+1} \\ 0, & \text{otherwise} \end{cases},\tag{22}$$

$$B^{i,j}(\mathbf{u}) = \frac{\mathbf{u} - \mathbf{u}^i}{\mathbf{u}^{i+j} - \mathbf{u}^i} B^{i,j-1}(\mathbf{u}) + \frac{\mathbf{u}^{i+j+1} - \mathbf{u}}{\mathbf{u}^{i+j+1} - \mathbf{u}^{i+1}} B^{i+1,j-1}(\mathbf{u}),\tag{23}$$

where **u***<sup>i</sup>* is the *i*th element of the knot vector for non-uniform open splines. Next, the numerical values of parameters, which are scalar, should be evaluated by using Equation (21) at time step *k* to be incorporated in Equation (20). The calculation of Equation (20) at the desired input signal is straightforward. However, the numerical computation of the Jacobian matrix **fx** *<sup>k</sup>*−<sup>1</sup> could be challenging when using a curve-fitting method. Each term of the Jacobian matrix can be approximated by using complex variables to reduce the truncation error [53,54] for very small increments. For instance, in the case of a multi-variable function, the Jacobian column of Equation (20) with respect to the *r*th term of the augmented state vector **x**ˆ *<sup>k</sup>*−<sup>1</sup> can be written in the partial derivative form as

$$\frac{\partial \mathbf{f}(\mathbf{\hat{x}}\_{k-1,1'}'\mathbf{\hat{x}}\_{k-1,2'}'\dots\mathbf{\hat{x}}\_{k-1,L}')}{\partial \mathbf{\hat{x}}\_{k-1,r}'} = \frac{Im(\mathbf{f}(\mathbf{\hat{x}}\_{k-1,1'}'\mathbf{\hat{x}}\_{k-1,2'}'\dots\mathbf{\hat{x}}\_{k-1,r}' + i\delta,\dots\mathbf{\hat{x}}\_{k-1,L}'))}{\delta} + O(\delta^2),\tag{24}$$

where *r* ∈ [1, *L*], and *iδ* represents a very small increment in the complex plane. *δ* = *Lε* is a real value. Epsilon *ε* is a very small number and depends on the machine's precision. The *r*th term of the state vector **x**ˆ *<sup>k</sup>*−1,*<sup>r</sup>* <sup>+</sup> *<sup>i</sup><sup>δ</sup>* is expanded by using the Taylor series [53,54].

The evaluation of *Im*(**f**(**x**ˆ<sup>+</sup> *<sup>k</sup>*−1,1, **<sup>x</sup>**ˆ<sup>+</sup> *<sup>k</sup>*−1,2, ..., **<sup>x</sup>**ˆ *<sup>k</sup>*−1,*<sup>r</sup>* <sup>+</sup> *<sup>i</sup>δ*, ..., **<sup>x</sup>**ˆ<sup>+</sup> *<sup>k</sup>*−1,*L*)) for the parameter vector **u***<sup>k</sup>* requires the evaluation of *C*(**uk**−**1**) as complex numbers. However, the knot vector cannot contain any complex numbers [48,49]. Therefore, a criterion |*t* − *t i <sup>k</sup>*−1<sup>|</sup> <sup>&</sup>lt; <sup>Ξ</sup> is introduced, where Ξ can be a small real number, such as 0.1, which implies that the knotpoint distance between *t i <sup>k</sup>*−<sup>1</sup> and the complex input argument *<sup>t</sup>* should be less than <sup>Ξ</sup>. Using this criterion enables the knot points to be evaluated with the curve-fitting method through complex input. With the Jacobian of the dynamic system **fx** *k*−1 , the covariance matrix **P**− *<sup>k</sup>* is approximated as [46]

$$\mathbf{P}\_{k}^{-} = \mathbf{f}\_{\mathbf{x}\_{k-1}^{\prime}} \mathbf{P}\_{k-1}^{+} \mathbf{f}\_{\mathbf{x}\_{k-1}^{\prime}}^{T} + \mathbf{Q}\_{k\prime} \tag{25}$$

where **Q***<sup>k</sup>* is the covariance matrix of the plant noise. In the update stage, the sensor measurements are taken into account to improve the estimated augmented state vector **x**ˆ− *k* . The innovation **Δ***<sup>k</sup>* is calculated as [46]

$$
\Delta\_k = \mathbf{o}\_k - \mathbf{h}(\mathbf{\hat{x}}\_{\mathbf{k}}^{\prime -}), \tag{26}
$$

where **o***<sup>k</sup>* are the sensor measurements at the *k* time step, and **h**(**x**− **<sup>k</sup>** ) is the sensor measurement function. With the Jacobian of the sensor measurement function **hx** , the innovation in the covariance matrix **S***<sup>k</sup>* and the Kalman gain **K***<sup>k</sup>* can be calculated as [46]

$$\begin{aligned} \mathbf{S}\_k &= \mathbf{h}\_{\mathbf{x}\_k'} \mathbf{P}\_k^{-} \mathbf{h}\_{\mathbf{x}\_k'}^T + \mathbf{R}\_k \\ \mathbf{K}\_k &= \mathbf{P}\_k^{-} \mathbf{h}\_{\mathbf{x}\_k'}^T \mathbf{S}\_k^{-1} \end{aligned} \qquad \begin{aligned} \mathbf{R}\_k &= \\ \end{aligned} $$

where **R***<sup>k</sup>* is the covariance matrix of the measurement noise. Finally, the augmented state vector **x**ˆ<sup>+</sup> *<sup>k</sup>* and covariance matrix **<sup>P</sup>**<sup>+</sup> *<sup>k</sup>* are updated at the time step *k*, which will be used for the next time step as [46]

$$\begin{aligned} \mathbf{\hat{x}}\_{k}^{\prime +} &= \mathbf{\hat{x}}\_{k}^{\prime -} + \mathbf{K}\_{k} \mathbf{\hat{A}}\_{k} \\ \mathbf{P}\_{k}^{+} &= (\mathbf{I}\_{L} - \mathbf{K}\_{k} \mathbf{h}\_{\mathbf{x}\_{k}^{\prime}}) \mathbf{P}\_{k}^{-} \end{aligned} \qquad \begin{aligned} \mathbf{\hat{x}}\_{k}^{\prime} &= \mathbf{\hat{x}}\_{k}^{\prime +} \\ \mathbf{P}\_{k}^{-} &= \mathbf{\hat{x}}\_{k}^{\prime +} \end{aligned} \tag{28}$$

where **I***<sup>L</sup>* is the identity matrix of dimension *L*.

#### Covariance Matrices of Process and Measurement Noises

It is well known that when applying Kalman filters, the tuning of the filter parameters is crucial, especially the covariance matrices of the plant and measurement noises. Furthermore, it was established in [30,31] that in dealing with non-linear systems, the improper tuning/setting of these covariance matrices can make the algorithm unstable. In this study, the properties of measurement noise are precisely known because the measurements are built from a dynamic model (providing ground truth) with an addition of white Gaussian noise to replicate real sensors. Thus, the covariance matrix of measurement noise can be obtained. For instance, when using position and pressure sensors, the covariance matrix of the measurement noise, **R**, would then take the form [30,31,34]

$$\mathbf{R} = \begin{bmatrix} \left(\sigma\_s'\right)^2 \mathbf{I}\_n & \mathbf{0}\_{n \times n\_p} \\ \mathbf{0}\_{n\_p \times n} & \left(\sigma\_p'\right)^2 \mathbf{I}\_{n\_p} \end{bmatrix} \tag{29}$$

where *σ <sup>s</sup>* and *σ <sup>p</sup>* are the standard deviations of measurement noises at the position and pressure levels, respectively. In Equation (29), *n* is the number of actuator sensors and *np* is the number of pressure sensors. **I***n*, **I***np* , **0***n*×*np* , and **0***np*×*<sup>n</sup>* are the identity and zero matrices of corresponding orders, respectively. In the case of a multibody model along with positions and velocities as the state vector, the structure of the plant noise in the discrete-time frame was well established in [30,31] and can be written as

$$\mathbf{Q} = \begin{bmatrix} \frac{\sigma\_\mathbf{x}^2 \Delta t^3 \mathbf{I}\_{n\_f}}{3} & \frac{\sigma\_\mathbf{x}^2 \Delta t^2 \mathbf{I}\_{n\_f}}{2} \\ \frac{\sigma\_\mathbf{x}^2 \Delta t^2 \mathbf{I}\_{n\_f}}{2} & \sigma\_\mathbf{x}^2 \Delta t \mathbf{I}\_{n\_f} \end{bmatrix},\tag{30}$$

where Δ*t* is the size of the integration time step and *nf* is the number of degrees of freedom of the system. It should be noted that Equation (30) includes the variance at the acceleration level, *σ***x¨**, because of the acceleration errors arising from the inaccurate description of forces and mass distribution. Furthermore, the state vector in this study also includes the hydraulic pressures and the hydraulic parameters, along with the positions and velocities, and errors can occur at the pressure and parameter levels as well. Therefore, inspired by [34], the variance of hydraulic pressures, *σp*,*D*, and the variance of hydraulic parameters, *σh*p,D , can be directly incorporated as the diagonal elements in Equation (30). Accordingly, the structure of the covariance matrix of the plant noise, **Q**, in the parameter estimation can be written as

$$\mathbf{Q} = \begin{bmatrix} \frac{\sigma\_{\mathbf{x}}^{2} \Delta t^{2} \mathbf{I}\_{n\_{f}}}{3} & \frac{\sigma\_{\mathbf{x}}^{2} \Delta t^{2} \mathbf{I}\_{n\_{f}}}{2} & \mathbf{0}\_{n\_{f} \times n\_{p} + n\_{h\_{p}}} & \mathbf{0}\_{n\_{f} \times n\_{p} + n\_{p}} \\ \frac{\sigma\_{\mathbf{x}}^{2} \Delta t^{2} \mathbf{I}\_{n\_{f}}}{2} & \sigma\_{\mathbf{x}}^{2} \Delta t \mathbf{I}\_{n\_{f}} & \mathbf{0}\_{n\_{f} \times n\_{p} + n\_{h\_{p}}} & \mathbf{0}\_{n\_{f} \times n\_{p} + n\_{h\_{p}}} \\ \mathbf{0}\_{n\_{p} + n\_{h\_{p}} \times n\_{f}} & \mathbf{0}\_{n\_{p} + n\_{h\_{p}} \times n\_{f}} & \sigma\_{p,D}^{2} & \mathbf{0}\_{n\_{f} \times n\_{p} + n\_{h\_{p}}} \\ \mathbf{0}\_{n\_{p} + n\_{h\_{p}} \times n\_{f}} & \mathbf{0}\_{n\_{p} + n\_{h\_{p}} \times n\_{f}} & \mathbf{0}\_{n\_{p} + n\_{h\_{p}} \times n\_{f}} & \sigma\_{h\_{p,D}}^{2} \end{bmatrix}. \tag{31}$$

In this study, the integration errors are assumed to be negligible in comparison to the acceleration, pressure, and parameter errors.

#### **3. Case Example: Hydraulically Actuated System**

The parameter estimation methodology described in Section 2 is applied to estimate the characteristic curves at the *a*, *b*, *c*, and *d* ports of the 4/3 directional control valve shown in Figure 2. A four-bar mechanism actuated by a hydraulic circuit is presented in Figure 2. The dynamics of the mechanism are modelled using the semi-recursive and hydraulic lumped fluid theories, as described below.
