Estimating the Characteristic Curve of a Directional Control Valve in a Combined Multibody and Hydraulic System Using an Augmented Discrete Extended Kalman Filter

Abstract

The estimation of the parameters of a simulation model such that the model’s behaviour matches closely with reality can be a cumbersome task. This is due to the fact that a number of model parameters cannot be directly measured, and such parameters might change during the course of operation in a real system. Friction between different machine components is one example of these parameters. This can be due to a number of reasons, such as wear. Nevertheless, if one is able to accurately define all necessary parameters, essential information about the performance of the system machinery can be acquired. This information can be, in turn, utilised for product-specific tuning or predictive maintenance. To estimate parameters, the augmented discrete extended Kalman filter with a curve fitting method can be used, as demonstrated in this paper. In this study, the proposed estimation algorithm is applied to estimate the characteristic curves of a directional control valve in a four-bar mechanism actuated by a fluid power system. The mechanism is modelled by using the double-step semi-recursive multibody formulation, whereas the fluid power system under study is modelled by employing the lumped fluid theory. In practise, the characteristic curves of a directional control valve is described by three to six data control points of a third-order B-spline curve in the augmented discrete extended Kalman filter. The results demonstrate that the highly non-linear unknown characteristic curves can be estimated by using the proposed parameter estimation algorithm. It is also demonstrated that the root mean square error associated with the estimation of the characteristic curve is 0.08% with respect to the real model. In addition, all the errors in the estimated states and parameters of the system are within the 95% confidence interval. The estimation of the characteristic curve in a hydraulic valve can provide essential information for performance monitoring and maintenance applications.

1. Introduction

Multibody system dynamics (MBS) approaches enable the creation of the equations of motion that describe a mechanical system and relevant sub-components of complex mechanical systems [1,2]. The use of MBS leads to physics-based models that act as a single source of information [3] and represent the operation of an equivalent physical system in the real world [4]. The data generated by an MBS simulation model can be used to solve real-world problems throughout a product’s life cycle [5].

A physical system might have parameters that are difficult to estimate and that could accordingly create uncertainties in MBS models. In the real world, these parameters might be cumbersome or sometimes even impossible to measure directly due to economical limitations and sensor implementation difficulties. In addition, these parameters might change over time due to wear and other factors that come into play during operation. In some cases, parameters can only be interpreted from the manufacturer’s catalogues while not manifesting the current state of a product or differences in individual products due to manufacturing tolerances. Estimating these parameters can provide valuable information about the state and working performance of a product [6,7]. Manufacturers can use this information for condition monitoring [8,9], predictive maintenance [10,11,12], and real-time simulations for digital-twin applications [13,14].

In general, parameter estimation is a discipline that provides the essential tools for the estimation of parameters appearing in the modelling of a system [15]. The most common techniques for parameter estimation are weighted least squares [16,17], Kalman filtering [18,19], orthogonal least squares [20], robust techniques that include clustering [10], and regression diagnostics [21]. Among these algorithms, Kalman filters for parameter estimation have been utilised in a wide variety of engineering studies, ranging from control [22] and mechatronics [23] to heat transfer [24,25], fluid mechanics [26,27], turbulence [28], and others.

In the MBS field, several types of Kalman filter algorithms have been used to estimate system states based on the multibody equations of motion [29,30,31,32,33,34]. In state estimation, the independent coordinate method was introduced by using the independent positions and velocities of the multibody model as the states of the Kalman filter [30]. Using the independent coordinate method, the MBS formulation offers a general approach for estimating the system coordinates in terms of independent states for open- and closed-loop systems [30,31,34]. Less attention has been given to parameter estimation in MBS systems [35]. This is due to the complexity of the problem. As in many cases, the parameters are not constant and have to be estimated from the measurable variables of the dynamic system. In an-MBS related study, vehicle suspended mass and road friction were estimated in a dual-estimation application by using the extended Kalman filter (EKF) and the unscented Kalman filter (UKF) [36]. The generalised polynomial chaos (gPC) theory was first implemented in the framework of MBS in 2006 to quantify the parametric and external uncertainties [37,38]. However, in [37,38], only constant parameters were estimated.

Contrary to this, in practical systems, the parameters are a function of several system variables and may follow very complicated and unknown non-linear variations during the working cycles [39,40]. For instance, in the case of a hydraulically actuated mobile working machine, the characteristic curve of a hydraulic valve can play a significant role in terms of machine performance [41,42]. The characteristic curve of a hydraulic valve can be expressed as a function of the spool position and the semi-empiric flow rate coefficient [41,42]. The semi-empiric flow rate coefficient relates the discharge coefficient, pressure losses, and flow characteristics that demonstrate the dynamic characteristics of a hydraulic valve [43,44,45]. Accordingly, the characteristic curve of a hydraulic valve can be used in the condition monitoring and predictive maintenance of hydraulically driven systems [42]. However, in an operating hydraulic system, only the minimum and maximum points on the characteristic curve can be determined from the manufacturer’s catalogues with a high level of certainty [41,42]. The characteristic curve of a hydraulic valve remains unclear in a working cycle and varies from one hydraulic valve to another due to manufacturing tolerances and possible wear [41,42]. Applying parameter estimation theories [46,47] in combination with MBS equations of motion can enable the estimation of the characteristic curve of a hydraulically driven physical system in operation by using a limited amount of information.

Generally, unknown parameters are treated as constants in the dynamic equations of motion. The estimation of non-linear parameters typically requires an accurate description of the first derivatives of the corresponding parameters. However, in the real world, the first derivatives of parameters are unclear. The first derivative of a characteristic curve in a hydraulic valve is an example. In the case of a characteristic curve, a vector of data points can be constructed using random points between the minimum and maximum values provided in the manufacturer’s catalogues. Through a parameter vector, the characteristic curve of a hydraulic valve in the real world can be estimated by combining parameter estimation algorithms [46,47] and curve-fitting methods [48,49,50,51,52]. Considering parameter estimation constraints, this study proposes the estimation of parameters by combining the augmented discrete extended Kalman filter (ADEKF) with curve-fitting methods.

The objective of this study is to propose a parameter estimation algorithm by combining the ADEKF algorithm with a curve-fitting method in an application for estimating linear and non-linear parameters. To this end, parameters are introduced as vectors in the augmented state vector. Due to the accuracy of the finite difference schemes in the complex plane, as demonstrated in [53,54], an approach to computing the Jacobian of a non-linear system of ordinary differential equations (ODEs) through complex variables in the framework of a parameter estimation algorithm is proposed. Based on the parameter estimation algorithm, the structures of covariance matrices of plant and measurement noises are introduced. The parameter estimation algorithm is applied to estimate the characteristic curve of a directional control valve in a hydraulically driven four-bar mechanism. As reported in [55], the double-step formulation has advantages over Index-3 Augmented Lagrangian formulation due to the use of a coordinate partitioning method [56]. Therefore, the double-step semi-recursive formulation is used to model the four-bar mechanism with relative coordinates. A fluid power system, in turn, is modelled by using the lumped fluid theory. This algorithm is verified by estimating the characteristic curves of the directional control valve using three, four, five, and six vector data control points in the mechanism. The implementation of the parameter estimation algorithm is explained by using MBS simulation models that represent the real model, estimation model, and simulation model. The estimation model considers the actuator position, pump pressure, and the pressure on the piston side as sensor measurements to account for the system responses. Applying the proposed parameter estimation methodology in MBS systems can enable the estimation of parameters of any complex system in a real-world system.

The rest of this paper is organised as follows. In Section 2, the parameter estimation methodology is described. Section 2 details further into the double-step semi-recursive MBS formulation, lumped fluid theory, monolithic approach, the ADEKF with a curve-fitting method, and structure of covariance matrices of plant and measurement noises. The parameter estimation methodology is applied to the case example presented in Section 3. Section 4 demonstrates the results of the parameter estimation algorithm for the case example. Finally, conclusions about parameter estimation are provided in Section 5.

2. Parameter Estimation Methodology

Figure 1 depicts a methodology that can be used to estimate the parameters of a dynamic system by using a simulation model. In this model, an initial covariance matrix ${\mathbf{P}}_{k-1}^{+}\in {\mathbb{R}}^{L\times L}$ and an augmented state vector ${\widehat{\mathbf{x}}}_{k-1}^{\prime +}={\left[\begin{array}{cc}{\mathbf{x}}_{k-1}^T& {\mathbf{y}}_{k-1}^T\end{array}\right]}^T$ at the time step $k-1$ are introduced. Here, L is the dimension of the augmented state vector, and $\mathbb{R}$ denotes the set of real numbers. $\mathbf{x}\in {\mathbb{R}}^{L-n_{h_p}}$ and $\mathbf{y}\in {\mathbb{R}}^{n_{h_p}}$ represent the states and parameters of the system, respectively. Here, $n_{h_p}$ is the number of hydraulic parameters.

In the real world, the sensors shown in the Figure 1 can be replaced by sensor measurements obtained from a physical system, such as a forklift, a tractor, etc. [4,57]. To account for the system response, the sensor measurement vector $\mathbf{o}$ includes the minimum number of measurements required by the ADEKF algorithm to estimate the states and parameters of a real system. In Figure 1, $\mathbf{h}$ corresponds to the sensor measurement function. Note that the parameters should not be included in the measurement vector, i.e., $\mathbf{y}\notin \mathbf{o}$. The parameter estimation algorithm estimates the augmented state vector ${\widehat{\mathbf{x}}}_k^{\prime +}$ and covariance matrix ${\mathbf{P}}_k^{+}$ from the minimum information of the real system at time step $k-1$ in the simulation model.

2.1. Multibody Dynamic Formulations

The parameter estimation methodology described in Section 2 is applied to the simulation of a hydraulically driven mechanism. In this study, the hydraulically driven mechanism is modelled using a double-step semi-recursive MBS formulation and the lumped fluid theory. The coupled multibody and hydraulic dynamics are integrated by using a single-step implicit trapezoidal integration scheme in a monolithic coupling approach. As demonstrated in [55], the double-step semi-recursive formulation uses a coordinate partitioning method [58,59,60] to express the hydraulically driven mechanism in terms of independent coordinates. As a result, the double-step semi-recursive formulation presents an appropriate multibody simulation approach for state and parameter estimation applications.

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 ${\mathbf{Z}}_i={\left[\begin{array}{cc}{\dot{\mathbf{r}}}_i^T& {\mathit{\omega }}_i^T\end{array}\right]}^T$ and six Cartesian accelerations as ${\dot{\mathbf{Z}}}_i={\left[\begin{array}{cc}{\ddot{\mathbf{r}}}_i^T& {\dot{\mathit{\omega }}}_i^T\end{array}\right]}^T$ for a complete description [61,62]. Here, ${\dot{\mathbf{r}}}_i$, ${\ddot{\mathbf{r}}}_i$, ${\mathit{\omega }}_i$ and ${\dot{\mathit{\omega }}}_i$ are velocities, accelerations, angular velocities, and angular accelerations of the body, respectively. In the relative coordinate system, the position of $n_b$ bodies in a system can be described by using joint coordinates as $\mathbf{z}={\left[\begin{array}{cccc}{z}_1& z_2& \dots & z_{n_b}\end{array}\right]}^T$ [61,62]. The absolute velocity $\mathbf{Z}$ and the absolute acceleration $\dot{\mathbf{Z}}$ of the system bodies can be mapped in terms of the relative joint velocity vector $\dot{\mathbf{z}}$ and the relative joint acceleration vector $\ddot{\mathbf{z}}$ by using the velocity transformation matrix as follows [61,62]:

$$\left.\begin{array}{cc}\hfill \mathbf{Z}& =\mathbf{T}{\mathbf{R}}_d\dot{\mathbf{z}}\hfill \\ \hfill \dot{\mathbf{Z}}& =\mathbf{T}{\mathbf{R}}_d\ddot{\mathbf{z}}+\mathbf{T}{\dot{\mathbf{R}}}_d\dot{\mathbf{z}}\hfill \end{array}\right\},$$

(1)

where $\mathbf{T}\in {\mathbb{R}}^{6n_b\times 6n_b}$ is the path matrix that demonstrates the topology of the system, and ${\mathbf{R}}_d\in {\mathbb{R}}^{6n_b\times n_b}$ is a block diagonal matrix. The path matrix $\mathbf{T}$ is a lower triangular matrix and contains entries of $6\times 6$ $\left({\mathbf{I}}_6\right)$ unit matrices representing bodies between the body under observation and the root of the system [61]. In Equation (1), the block diagonal matrix ${\mathbf{R}}_d$ and the product ${\dot{\mathbf{R}}}_d\dot{\mathbf{z}}$ can be computed with the joint-dependent element of the velocity transformation matrix ${\mathbf{b}}_i\in {\mathbb{R}}^{6\times 1}$ and the joint-dependent element of the acceleration transformation vector ${\mathbf{d}}_i\in {\mathbb{R}}^{6\times 1}$, respectively [61,62]. The semi-recursive formulation is described hereafter, but the interested reader is referred to [61,62] for further details of $\mathbf{T}$, ${\mathbf{b}}_i$, and ${\mathbf{d}}_i$. The composite mass matrix ${\overline{\mathbf{M}}}_i\in {\mathbb{R}}^{6\times 6}$ and the composite force vector ${\overline{\mathbf{Q}}}_i\in {\mathbb{R}}^{6\times 1}$ of the ith body can be described using absolute coordinates as [61,62]:

$${\overline{\mathbf{M}}}_i=\left[\begin{array}{cc}{m}_i{\mathbf{I}}_3& -m_i{\tilde{\mathbf{g}}}_i\\ m_i{\tilde{\mathbf{g}}}_i& {\mathbf{J}}_i-m_i{\tilde{\mathbf{g}}}_i{\tilde{\mathbf{g}}}_i\end{array}\right],$$

(2)

$${\overline{\mathbf{Q}}}_i=\left[\begin{array}{c}{\overline{\mathbf{f}}}_i-{\tilde{\mathit{\omega }}}_i\left({\tilde{\mathit{\omega }}}_i{m}_i{\mathbf{g}}_i\right)\\ {\mathbf{n}}_i-{\tilde{\mathit{\omega }}}_i{\mathbf{J}}_i{\mathit{\omega }}_i+{\tilde{\mathbf{g}}}_i({\mathbf{f}}_i-{\tilde{\mathit{\omega }}}_i\left({\tilde{\mathit{\omega }}}_i{m}_i{\mathbf{g}}_i\right))\end{array}\right],$$

(3)

where $m_i$ is the mass of the i body, ${\overline{\mathbf{f}}}_i\in {\mathbb{R}}^{3\times 1}$ is a vector of external forces, ${\tilde{\mathit{\omega }}}_i$ is the skew-symmetric matrix of the angular velocity vector, ${\mathbf{n}}_i\in {\mathbb{R}}^{3\times 1}$ is the vector of external moments, ${\mathbf{I}}_3$ is a $3\times 3$ unit matrix, and ${\tilde{\mathbf{g}}}_i\in {\mathbb{R}}^{3\times 3}$ 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), ${\mathbf{J}}_i$ is the inertia tensor of the ith 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}^T{\mathbf{T}}^T\overline{\mathbf{M}}\mathbf{T}{\mathbf{R}}_d\ddot{\mathbf{z}}={{\mathbf{R}}_d}^T{\mathbf{T}}^T(\overline{\mathbf{Q}}-\overline{\mathbf{M}}\mathbf{T}{\dot{\mathbf{R}}}_d\dot{\mathbf{z}}),$$

(4)

where $\overline{\mathbf{M}}\in {\mathbb{R}}^{6n_b\times 6n_b}$ is the block diagonal matrix consisting of the composite mass matrices of the bodies. The force vector $\overline{\mathbf{Q}}\in {\mathbb{R}}^{6n_b\times 1}$ 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 $\mathsf{\Phi }\left(\mathbf{z}\right) = \mathbf{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 ${\mathsf{\Phi }}_{\mathbf{z}}$ [62,63,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 $\dot{\mathbf{z}}$ can be described using the coordinate partitioning method [59]:

$$\left[\begin{array}{c}{\dot{\mathbf{z}}}^{\mathrm{d}}\\ {\dot{\mathbf{z}}}^{\mathrm{i}}\end{array}\right]=\left[\begin{array}{c}-{\left({\mathsf{\Phi }}_{\mathbf{z}}^{\mathrm{d}}\right)}^{-1}{\mathsf{\Phi }}_{\mathbf{z}}^{\mathrm{i}}\\ \mathbf{I}\end{array}\right]{\dot{\mathbf{z}}}^{\mathrm{i}}\equiv {\mathbf{R}}_{\mathrm{z}}{\dot{\mathbf{z}}}^{\mathrm{i}},$$

(5)

where ${\dot{\mathbf{z}}}^{\mathrm{d}}\in {\mathbb{R}}^m$ are the relative joint-dependent velocities, ${\dot{\mathbf{z}}}^{\mathrm{i}}\in {\mathbb{R}}^{n_f}$ are the relative joint-independent velocities, ${\mathbf{R}}_{\mathrm{z}}\in {\mathbb{R}}^{n_b\times n_f}$ is the velocity transformation matrix, ${\mathsf{\Phi }}_{\mathbf{z}}^{\mathrm{d}}\in {\mathbb{R}}^{m\times m}$, and ${\mathsf{\Phi }}_{\mathbf{z}}^{\mathrm{i}}\in {\mathbb{R}}^{m\times n_f}$ 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 ${\mathsf{\Phi }}_{\mathbf{z}}^{\mathrm{d}}$ can be obtained [55,59]. The relative joint acceleration vector can be expressed by differentiating Equation (5) [59]:

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

(6)

where ${\ddot{\mathbf{z}}}^{\mathrm{i}}$ are the relative joint-independent accelerations, and ${\dot{\mathbf{R}}}_{\mathrm{z}}$ 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}}_{\mathrm{z}}^T{\mathbf{R}}_d^T{\mathbf{T}}^T\overline{\mathbf{M}}\mathbf{T}{\mathbf{R}}_d{\mathbf{R}}_{\mathrm{z}}{\ddot{\mathbf{z}}}^{\mathrm{i}}={\mathbf{R}}_{\mathrm{z}}^T{\mathbf{R}}_d^T\left({\mathbf{T}}^T\overline{\mathbf{Q}}-{\mathbf{T}}^T\overline{\mathbf{M}}\mathbf{D}\right),$$

(7)

where $\mathbf{D}=\mathbf{T}{\mathbf{R}}_d\left[\begin{array}{c}-{\left({\mathsf{\Phi }}_{\mathbf{z}}^{\mathrm{d}}\right)}^{-1}\left({\dot{\mathsf{\Phi }}}_{\mathbf{z}}\dot{\mathbf{z}}\right)\\ \mathbf{0}\end{array}\right]+\mathbf{T}{\dot{\mathbf{R}}}_d\dot{\mathbf{z}}$ represent the absolute accelerations when the vector $\ddot{\mathbf{z}}$ is zero in Equation (6). Equation (7) can be further simplified using the accumulated mass matrix ${\mathbf{M}}^{\Sigma }={\mathbf{R}}_{\mathrm{z}}^T{\mathbf{R}}_d^T{\mathbf{T}}^T\overline{\mathbf{M}}\mathbf{T}{\mathbf{R}}_d{\mathbf{R}}_{\mathrm{z}}$ and the accumulated force matrix ${\mathbf{Q}}^{\Sigma }={\mathbf{R}}_{\mathrm{z}}^T{\mathbf{R}}_d^T\left({\mathbf{T}}^T\overline{\mathbf{Q}}-{\mathbf{T}}^T\overline{\mathbf{M}}\mathbf{D}\right)$.

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 $V_h$, the differential pressure ${\dot{p}}_h$ can be computed [41,65] as

$${\dot{p}}_h=\frac{k_p+p_h{k}_0}{V_h}\sum _{h=1}^{n_f}{Q}_h,$$

(8)

where $Q_h$ is the sum of incoming and outgoing volume flow rates, $k_0$ is the flow gain, $k_p$ is the pressure flow coefficient, and $n_f$ is the total amount of volume flows. By employing a semi-empirical method, the volume flow rate $Q_R$ through a throttle valve can be described as [66]

$$Q_R=C_R\mathrm{sgn}\left(\Delta p\right)\sqrt{\mid \Delta p\mid },$$

(9)

where $\Delta p$ is the pressure difference and $C_R=C_d{A}_R\sqrt{\frac{2}{\rho }}$ is the semi-empirical flow rate coefficient of the throttle valve. Here, $C_d$ is the flow discharge coefficient and $\rho$ is the fluid density. In Equation (9), $A_R$ is the cross-sectional area of the pressure relief valve. Similarly, the volume flow rate $Q_d$ through a directional control valve can be computed as [67,68]

$$Q_d=C_{v}U\mathrm{sgn}\left(\Delta p\right)\sqrt{\mid \Delta p\mid },$$

(10)

where $C_v$ 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{U}=\frac{U_{ref}-U}{\tau },$$

(11)

where $U_{ref}$ is the reference voltage signal, and $\tau$ is the time constant. The incoming flow rate $Q_{in}$ and outgoing flow rate $Q_{out}$ in the hydraulic cylinder can be described as

$$Q_{in}=\dot{s}{A}_1,Q_{out}=\dot{s}{A}_2,$$

(12)

where $\dot{s}$ is the piston velocity, and $A_1$ and $A_2$ are the areas on the piston and piston-rod side of the cylinder, respectively. The force $F_h$ produced by the cylinder can be written as

$$F_h=p_1{A}_1-p_2{A}_2-F_{\mu },$$

(13)

where $p_1$ and $p_2$ are, respectively, the pressure on the piston and piston-rod side, which can be calculated by using Equation (8). $F_{\mu }$ 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{\Vert \dot{\mathbf{s}}\Vert }{{\mathrm{v}}_s}\right)+(F_s-F_c)\frac{\frac{\Vert \dot{\mathbf{s}}\Vert }{{\mathrm{v}}_s}}{{\left(\begin{array}{c}\frac{1}{4}{\left(\frac{\Vert \dot{\mathbf{s}}\Vert }{{\mathrm{v}}_s}\right)}^2+\frac{3}{4}\end{array}\right)}^2}\right)\mathrm{sgn}\left(\dot{\mathbf{s}}\right)+{\sigma }_2\dot{\mathbf{s}}tanh\left(4\right),$$

(14)

where $F_c$ is the Coulomb friction, ${\mathrm{v}}_s$ is the Stribeck velocity, $F_s$ is the static friction, ${\sigma }_2$ is the coefficient of viscous friction, and $\dot{\mathbf{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:

$$\left.\begin{array}{cc}\hfill {\mathbf{M}}^{\Sigma }\left(\mathbf{z}\right){\ddot{\mathbf{z}}}^{\mathrm{i}}& ={\mathbf{Q}}^{\Sigma }(\mathbf{z},\dot{\mathbf{z}},\mathbf{p})\hfill \\ \hfill \dot{\mathbf{p}}& =\mathbf{v}(\mathbf{z},\dot{\mathbf{z}},\mathbf{p},\mathbf{y},U)\hfill \end{array}\right\},$$

(15)

where $\mathbf{p}$ is the pressure vector, and $\mathbf{y}$ is the vector of hydraulic parameters. Equation (15) is a set of non-linear equations that can be represented as $\mathbf{f}(\overline{\mathbf{x}},U)=\mathbf{0}$. Here, the vector $\overline{\mathbf{x}}={\left[\begin{array}{cccc}{\mathbf{z}}^T& {\dot{\mathbf{z}}}^T& {\mathbf{p}}^T& {\mathbf{y}}^T\end{array}\right]}^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,72,73]. In this integration scheme, the relative joint-independent positions and the pressures are initially predicted as ${\mathbf{z}}_{k+1}^{\mathrm{i}}={\mathbf{z}}_k^{\mathrm{i}}+{\dot{\mathbf{z}}}_k^{\mathrm{i}}\Delta k+\frac{1}{2}{\ddot{\mathbf{z}}}_k^{\mathrm{i}}\Delta k^2$ and ${\mathbf{p}}_{k+1}={\mathbf{p}}_k+{\dot{\mathbf{p}}}_k\Delta k$, respectively [73]. The derivatives of ${\mathbf{z}}_{k+1}^{\mathrm{i}}$ and ${\mathbf{p}}_{k+1}$ can be predicted as

$$\left.\begin{array}{cc}\hfill {\dot{\mathbf{z}}}_{k+1}^{\mathrm{i}}& =\frac{2}{\Delta k}{\mathbf{z}}_{k+1}^{\mathrm{i}}+{\check{\dot{\mathbf{z}}}}_k^{\mathrm{i}}\hfill \\ \hfill {\ddot{\mathbf{z}}}_{k+1}^{\mathrm{i}}& =\frac{4}{\Delta k^2}{\mathbf{z}}_{k+1}^{\mathrm{i}}+{\check{\ddot{\mathbf{z}}}}_k^{\mathrm{i}}\hfill \\ \hfill {\dot{\mathbf{p}}}_{k+1}& =\frac{2}{\Delta k}{\mathbf{p}}_{k+1}+{\check{\dot{\mathbf{p}}}}_k\hfill \end{array}\right\},$$

(16)

where ${\check{\dot{\mathbf{z}}}}_k^{\mathrm{i}}=-(\frac{2}{\Delta k}{\mathbf{z}}_k^{\mathrm{i}}+{\dot{\mathbf{z}}}_k^{\mathrm{i}})$, ${\check{\ddot{\mathbf{z}}}}_k^{\mathrm{i}}=-(\frac{4}{\Delta k^2}{\mathbf{z}}_k^{\mathrm{i}}+\frac{4}{\Delta k}{\dot{\mathbf{z}}}_k^{\mathrm{i}}+{\ddot{\mathbf{z}}}_k^{\mathrm{i}})$ and ${\check{\dot{\mathbf{p}}}}_k=-(\frac{2}{\Delta k}{\mathbf{p}}_k+{\dot{\mathbf{p}}}_k)$. Note that the relative joint-dependent positions ${\mathbf{z}}_{k+1}^{\mathrm{d}}$ are obtained from ${\mathbf{z}}_{k+1}^{\mathrm{i}}$ and the previous step ${\mathbf{z}}_k^{\mathrm{d}}$ by solving the position problem $\mathsf{\Phi }\left(\mathbf{z}\right) = 0$ [59,60,61]. The non-linear constraint equations are solved iteratively with the Newton–Raphson method [59,60,61]. The derivatives of the relative joint-dependent positions ${\mathbf{z}}_{k+1}^{\mathrm{d}}$ 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 $\overline{\mathbf{F}}\left({\chi }_{k+1}\right)=\mathbf{0}$ at the time step $k+1$ as

$$\left.\begin{array}{cc}\hfill {\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 }{\check{\ddot{\mathbf{z}}}}_k^{\mathrm{i}}& =\mathbf{0}\hfill \\ \hfill \frac{\Delta k}{2}{\mathbf{p}}_{k+1}-\frac{\Delta k^2}{4}{\mathbf{v}}_{k+1}+\frac{\Delta k^2}{4}{\check{\dot{\mathbf{p}}}}_k& =\mathbf{0}\hfill \end{array}\right\},$$

(17)

where ${\chi }_{k+1}={\left[\begin{array}{cc}{\left({\mathbf{z}}^{\mathrm{i}}\right)}_{k+1}^T& {\mathbf{p}}_{k+1}^T\end{array}\right]}^T$ is unknown. The Newton–Raphson method is employed on the non-linear Equation (17) to iteratively compute the unknown variables [73,74]:

$${\left[\begin{array}{c}\frac{\partial \overline{\mathbf{F}}\left(\chi \right)}{\partial \chi }\end{array}\right]}_{k+1}^{\left(h\right)}\Delta {\chi }_{k+1}^{\left(h\right)}=-{\left[\begin{array}{c}\overline{\mathbf{F}}\left(\chi \right)\end{array}\right]}_{k+1}^{\left(h\right)},$$

(18)

where $\Vert \Delta {\mathbf{z}}_{k+1}^{\mathrm{i}}\Vert <1\times {10}^{-10}\mathrm{rad}$ and $\Vert \Delta {\mathbf{p}}_{k+1}\Vert <1\times {10}^{-2}\mathrm{Pa}$ are the error tolerances in the relative joint independent positions and pressures provided in the Newton–Raphson method. In Equation (18), ${\left[\begin{array}{c}\overline{\mathbf{F}}\left(\chi \right)\end{array}\right]}_{k+1}^{\left(h\right)}$ is the residual vector, which can be computed as

$${\left[\begin{array}{c}\overline{\mathbf{F}}\left(\chi \right)\end{array}\right]}_{k+1}^{\left(h\right)}=\frac{\Delta k^2}{4}{\left[\begin{array}{c}{\mathbf{M}}^{\Sigma }{\ddot{\mathbf{z}}}^{\mathrm{i}}-{\mathbf{Q}}^{\Sigma }\\ \dot{\mathbf{p}}-\mathbf{v}\end{array}\right]}_{k+1}^{\left(h\right)}.$$

(19)

In Equation (18), ${\left[\begin{array}{c}\frac{\partial \overline{\mathbf{F}}\left(\chi \right)}{\partial \chi }\end{array}\right]}_{k+1}^{\left(h\right)}$ is the tangent matrix, which can be numerically approximated at a point ${\chi }_0$ by using a forward finite differences, as demonstrated in the literature [72,75]. Now, by computing $\Delta {\chi }_{k+1}^{(h+1)}$ from Equation (18), the next iteration ${\chi }_{k+1}^{(h+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 ${\mathbf{x}}^{\prime }={\left[\begin{array}{cccc}{\left({\mathbf{z}}^{\mathrm{i}}\right)}^T& {\left({\dot{\mathbf{z}}}^{\mathrm{i}}\right)}^T& {\mathbf{p}}^T& {\mathbf{y}}^T\end{array}\right]}^T$. At this step, ${\widehat{\mathbf{x}}}_k^{\prime -}$ is calculated in the time integration of a dynamic model described as [46]

$${\widehat{\mathbf{x}}}_k^{\prime -}=\mathbf{f}({\widehat{\mathbf{x}}}_{k-1}^{\prime +},U_k),$$

(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 $\mathbf{u}$ for non-uniform open splines [48,49] at the current time step:

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

(21)

where n is the number of control points, d is the degree, $B^{i,d}\left(\mathbf{u}\right)$ are the dth order of B-spline basis functions, and ${\mathbf{N}}^i$ is the control point vector. The control point vector can be expressed in terms of the system parameters $\mathbf{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 $\mathbf{N}=\left[\begin{array}{ccccc}{U}_{\min }& U_1& \dots & U_{\mathrm{n}}& U_{\max }\\ C_{v_{\min }}& C_{v_1}& \dots & C_{v_{\mathrm{n}}}& C_{v_{\max }}\end{array}\right]$. Here, $U_{\min },U_1$, and $U_{\mathrm{n}}$ represent spool positions, and $C_{v_{\min }},C_{v_1}$, and $C_{v_{\max }}$ are the semi-empiric flow rate coefficients of a hydraulic valve. $B^{i,d}\left(\mathbf{u}\right)$ can be defined by using the Cox–de Boor recursion formula [48,49]:

$$B^{i,0}\left(\mathbf{u}\right)=\left\{\begin{array}{cc}1\hfill & {\mathbf{u}}^i\le \mathbf{u}<{\mathbf{u}}^{i+1}\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.,$$

(22)

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

(23)

where ${\mathbf{u}}^i$ is the ith 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 ${\mathbf{f}}_{{\mathbf{x}}_{k-1}^{\prime }}$ 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 rth term of the augmented state vector ${{\widehat{\mathbf{x}}}^{\prime }}_{k-1}$ can be written in the partial derivative form as

$$\frac{\partial \mathbf{f}({\widehat{\mathbf{x}}}_{k-1,1}^{\prime },{\widehat{\mathbf{x}}}_{k-1,2}^{\prime },\dots ,{\widehat{\mathbf{x}}}_{k-1,L}^{\prime })}{\partial {\widehat{\mathbf{x}}}_{k-1,r}^{\prime }}=\frac{Im\left(\mathbf{f}({\widehat{\mathbf{x}}}_{k-1,1}^{\prime +},{\widehat{\mathbf{x}}}_{k-1,2}^{\prime +},\dots ,{\widehat{\mathbf{x}}}_{k-1,r}^{\prime }+i\delta ,\dots ,{\widehat{\mathbf{x}}}_{k-1,L}^{\prime +})\right)}{\delta }+O\left({\delta }^2\right),$$

(24)

where $r\in \left[1,L\right]$, and $i\delta$ represents a very small increment in the complex plane. $\delta =L\epsilon$ is a real value. Epsilon $\epsilon$ is a very small number and depends on the machine’s precision. The rth term of the state vector ${\widehat{\mathbf{x}}}_{k-1,r}^{\prime }+i\delta$ is expanded by using the Taylor series [53,54].

The evaluation of $Im\left(\mathbf{f}({\widehat{\mathbf{x}}}_{k-1,1}^{\prime +},{\widehat{\mathbf{x}}}_{k-1,2}^{\prime +},\dots ,{\widehat{\mathbf{x}}}_{k-1,r}^{\prime }+i\delta ,\dots ,{\widehat{\mathbf{x}}}_{k-1,L}^{\prime +})\right)$ for the parameter vector ${\mathbf{u}}_k$ requires the evaluation of $C\left({\mathbf{u}}_{\mathbf{k}-\mathbf{1}}\right)$ as complex numbers. However, the knot vector cannot contain any complex numbers [48,49]. Therefore, a criterion $|t-t_{k-1}^i|<\Xi$ is introduced, where $\Xi$ can be a small real number, such as $0.1$, which implies that the knot-point distance between $t_{k-1}^i$ and the complex input argument t should be less than $\Xi$. 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 ${\mathbf{f}}_{{\mathbf{x}}_{k-1}^{\prime }}$, the covariance matrix ${\mathbf{P}}_k^{-}$ 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,$$

(25)

where ${\mathbf{Q}}_k$ 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 ${\widehat{\mathbf{x}}}_k^{\prime -}$. The innovation ${\Delta }_k$ is calculated as [46]

$${\Delta }_k={\mathbf{o}}_k-\mathbf{h}\left({\widehat{\mathbf{x}}}_{\mathbf{k}}^{\prime -}\right),$$

(26)

where ${\mathbf{o}}_k$ are the sensor measurements at the k time step, and $\mathbf{h}\left({\mathbf{x}}_{\mathbf{k}}^{\prime -}\right)$ is the sensor measurement function. With the Jacobian of the sensor measurement function ${\mathbf{h}}_{{\mathbf{x}}^{\prime }}$, the innovation in the covariance matrix ${\mathbf{S}}_k$ and the Kalman gain ${\mathbf{K}}_k$ can be calculated as [46]

$$\left.\begin{array}{cc}\hfill {\mathbf{S}}_k& ={\mathbf{h}}_{{\mathbf{x}}_k^{\prime }}{\mathbf{P}}_k^{-}{\mathbf{h}}_{{\mathbf{x}}_k^{\prime }}^T+{\mathbf{R}}_k\hfill \\ \hfill {\mathbf{K}}_k& ={\mathbf{P}}_k^{-}{\mathbf{h}}_{{\mathbf{x}}_k^{\prime }}^T{\mathbf{S}}_k^{-1}\hfill \end{array}\right\},$$

(27)

where ${\mathbf{R}}_k$ is the covariance matrix of the measurement noise. Finally, the augmented state vector ${\widehat{\mathbf{x}}}_k^{\prime +}$ and covariance matrix ${\mathbf{P}}_k^{+}$ are updated at the time step k, which will be used for the next time step as [46]

$$\left.\begin{array}{cc}\hfill {\widehat{\mathbf{x}}}_k^{\prime +}& ={\widehat{\mathbf{x}}}_k^{\prime -}+{\mathbf{K}}_k{\Delta }_k\hfill \\ \hfill {\mathbf{P}}_k^{+}& =({\mathbf{I}}_L-{\mathbf{K}}_k{\mathbf{h}}_{{\mathbf{x}}_k^{\prime }}){\mathbf{P}}_k^{-}\hfill \end{array}\right\},$$

(28)

where ${\mathbf{I}}_L$ 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, $\mathbf{R}$, would then take the form [30,31,34]

$$\mathbf{R}=\left[\begin{array}{cc}{\left(\begin{array}{c}{\sigma }_s^{\prime }\end{array}\right)}^2{\mathbf{I}}_n& {\mathbf{0}}_{n\times n_p}\\ {\mathbf{0}}_{n_p\times n}& {\left(\begin{array}{c}{\sigma }_p^{\prime }\end{array}\right)}^2{\mathbf{I}}_{n_p}\end{array}\right],$$

(29)

where ${\sigma }_s^{\prime }$ and ${\sigma }_p^{\prime }$ 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 $n_p$ is the number of pressure sensors. ${\mathbf{I}}_n$, ${\mathbf{I}}_{n_p}$, ${\mathbf{0}}_{n\times n_p}$, and ${\mathbf{0}}_{n_p\times n}$ 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}=\left[\begin{array}{cc}{\displaystyle \frac{{\sigma }_{\ddot{\mathbf{x}}}^2\Delta t^3{\mathbf{I}}_{n_f}}{3}}& {\displaystyle \frac{{\sigma }_{\ddot{\mathbf{x}}}^2\Delta t^2{\mathbf{I}}_{n_f}}{2}}\\ {\displaystyle \frac{{\sigma }_{\ddot{\mathbf{x}}}^2\Delta t^2{\mathbf{I}}_{n_f}}{2}}& {\sigma }_{\ddot{\mathbf{x}}}^2\Delta t{\mathbf{I}}_{n_f}& \end{array}\right],$$

(30)

where $\Delta t$ is the size of the integration time step and $n_f$ is the number of degrees of freedom of the system. It should be noted that Equation (30) includes the variance at the acceleration level, ${\sigma }_{\ddot{\mathbf{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, ${\sigma }_{p,D}$, and the variance of hydraulic parameters, ${\sigma }_{h_{\mathrm{p},\mathrm{D}}}$, can be directly incorporated as the diagonal elements in Equation (30). Accordingly, the structure of the covariance matrix of the plant noise, $\mathbf{Q}$, in the parameter estimation can be written as

$$\mathbf{Q}=\left[\begin{array}{cccc}{\displaystyle \frac{{\sigma }_{\ddot{\mathbf{x}}}^2\Delta t^3{\mathbf{I}}_{n_f}}{3}}& {\displaystyle \frac{{\sigma }_{\ddot{\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_{h_p}}\\ {\displaystyle \frac{{\sigma }_{\ddot{\mathbf{x}}}^2\Delta t^2{\mathbf{I}}_{n_f}}{2}}& {\sigma }_{\ddot{\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{array}\right].$$

(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.

3.1. Dynamic Model of the System

The bodies of the mechanism are assumed to be rectangular beams whose lengths are $L_1=2\mathrm{m}$, $L_2=8\mathrm{m}$, and $L_3=5\mathrm{m}$ and whose masses are $m_1=100\mathrm{k}\mathrm{g}$, $m_2=400\mathrm{k}\mathrm{g}$, and $m_3=250\mathrm{k}\mathrm{g}$, respectively. The position vector at point D is ${\mathbf{r}}_D={\left[\begin{array}{ccc}-\frac{L_1}{2}& 0& 0\end{array}\right]}^T$. The point G is located at the centre of mass of body 1. The double-step semi-recursive formulation described in Section 2.1.1 is used to model the four-bar mechanism.

The four-bar mechanism is actuated using the sinusoidal reference input signal, which is taken as $U_{ref}=10sin\left(0.4\pi k\right)$, where k is the simulation run time. The simulations are performed for 5 $\mathrm{s}$. The hydraulic circuit consists of a double-acting hydraulic cylinder, connecting hoses 1 and 2, a 4/3 directional control valve, a pressure relief valve, a connecting hose of volume $V_p$, a differential pump of pressure $p_p$, and a tank with a constant pressure source $p_T$.

The lumped fluid theory described in Section 2.1.2 is used to compute the pressures within the hydraulic circuit. In the application of the lumped fluid theory, the hydraulic circuit can be divided into three control volumes $V_p$, $V_1$, and $V_2$. The pressure derivatives ${\dot{p}}_p$, ${\dot{p}}_1$, and ${\dot{p}}_2$ through these volumes can be computed as

$$\left.\begin{array}{cc}\hfill {\dot{p}}_p& =\frac{k_p+p_p{k}_0}{V_p}(Q_p-Q_R-Q_{d_1})\hfill \\ \hfill {\dot{p}}_1& =\frac{k_p+p_1{k}_0}{V_1}(Q_{d_1}-A_1\dot{s})\hfill \\ \hfill {\dot{p}}_2& =\frac{k_p+p_2{k}_0}{V_2}(A_2\dot{s}-Q_{d_2})\hfill \end{array}\right\},$$

(32)

where $Q_{d_1}$ and $Q_{d_2}$ are the flow rates in the control volumes 1 and 2. In Equation (32), $Q_p$ and $Q_R$ are the pump flow rate and flow rate through the pressure relief valve, respectively. The flow rates $Q_R$, $Q_{d_1}$, and $Q_{d_2}$ can be computed by employing Equations (9) and (10), respectively. The constant hydraulic parameters are tabulated in

Table 1. Parameters of the hydraulic circuit.

Parameter	Symbol	Value
Pump flow rate	$Q_p$	0.001 ${\mathrm{m}}^3/\mathrm{s}$
Tank pressure	$p_T$	0.1 MPa
Volume of the hose p	$V_p$	3.42 $\times {10}^{-3}{\mathrm{m}}^3$
Volume of the hose 1	$V_{h_1}$	3.42 $\times {10}^{-1}{\mathrm{m}}^3$
Volume of the hose 2	$V_{h_2}$	3.42 $\times {10}^{-1}{\mathrm{m}}^3$
Oil density	$\rho$	869 $\mathrm{k}\mathrm{g}$/$\mathrm{m}$${}^3$
Hydraulic parameter	$k_p$	1600 MPa
Hydraulic parameter	$k_0$	0.5
Area of the piston	$A_1$	2 $\times {10}^{-3}{\mathrm{m}}^2$
Area of the piston-rod	$A_2$	1.8 $\times {10}^{-3}{\mathrm{m}}^2$
Length of the cylinder/piston	l	$\sqrt{3}$$\mathrm{m}$
Area of pressure relief valve	$A_r$	2.24 $\times {10}^{-12}{\mathrm{m}}^2$
Area of directional control valve	$A_d$	1.96 $\times {10}^{-6}{\mathrm{m}}^2$
Coulomb friction force	$F_c$	210 N
Static friction force	$F_s$	830 N
Stribeck velocity	$v_s$	1.25 $\times {10}^{-2}$ $\mathrm{m}$/$\mathrm{s}$
Coefficient of viscous friction	$\sigma$	330 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$
Discharge coefficient	$C_d$	0.5
Area of throttle	$A_R$	2.24 $\times {10}^{-12}{\mathrm{m}}^2$

. In Equation (32), $\dot{s}$ is the actuator velocity, which can be determined from the actuator position vector $\mathbf{s}$. Following Figure 2, the vector $\mathbf{s}$ can be calculated from the position vectors ${\mathbf{r}}_G$ and ${\mathbf{r}}_D$ as

$$\left.\begin{array}{cc}\hfill \mathbf{s}& ={\mathbf{r}}_G-{\mathbf{r}}_D\hfill \\ \hfill \dot{s}& =\frac{d\mid \mathbf{s}\mid }{\Delta t}=\dot{\mathbf{s}}·\frac{\mathbf{s}}{\mid \mathbf{s}\mid }={\dot{\mathbf{r}}}_G·\frac{\mathbf{s}}{\mid \mathbf{s}\mid }\hfill \end{array}\right\},$$

(33)

where ${\dot{\mathbf{r}}}_G$ is the velocity vector of point G. The control volumes $V_1$ and $V_2$ appearing in Equation (32) can be calculated as follows:

$$\left.\begin{array}{cc}\hfill V_1& =V_{h_1}+A_1{l}_1\hfill \\ \hfill V_2& =V_{h_2}+A_2{l}_2\hfill \end{array}\right\},$$

(34)

where $V_{h_1}$, $V_{h_2}$, and $V_p$ are the control volumes of the respective hoses, as described in Table 1. In Equation (34), $l_1$ and $l_2$ are the lengths of the piston side and the piston-rod side chambers, respectively. $l_1$ and $l_2$ can be calculated with the vector $\mathbf{s}$ as

$$\left.\begin{array}{cc}\hfill l_1& =l_{1_0}-\left|\mathbf{s}\right|+s_0\hfill \\ \hfill l_2& =l_{2_0}+\left|\mathbf{s}\right|-s_0\hfill \end{array}\right\},$$

(35)

where $l_{1_0}$ and $l_{2_0}$ are the initial piston side length and the initial piston-rod side length, respectively. $l_{1_0}$ and $l_{2_0}$ are computed from the length of cylinder l, which is given in Table 1.

Using the vector $\mathbf{s}$, the hydraulic force ${\mathbf{F}}_h$ produced by the double-acting cylinder can be calculated as

$${\mathbf{F}}_h={\left[\begin{array}{ccc}\frac{s_X}{\mid \mathbf{s}\mid }{F}_h& \frac{s_Y}{\mid \mathbf{s}\mid }{F}_h& \frac{s_Z}{\mid \mathbf{h}\mid }{F}_h\end{array}\right]}^T,$$

(36)

where $F_h$ is computed from Equation (13). The hydraulic force vector ${\mathbf{F}}_h$ is combined with the external force vector ${\mathbf{f}}_i$ to calculate $\overline{\mathbf{Q}}$ in Equation (3). The resultant equations of motion (15) are formulated for the hydraulically driven four-bar mechanism. Equations (15) are solved by using an implicit single-step trapezoidal integration scheme in a monolithic approach, which was described in Section 2.1.3.

3.1.1. Real and Estimation Models

In this study, three dynamic versions of the mechanism are used to demonstrate the implementation of the parameter estimation algorithm. One of the models is the real model. The sensor measurements $\mathbf{o}$ are taken from the real model. The modelling errors are introduced in the force model of the estimation model with respect to the real model. The properties of the estimation model and the simulation model are the same. In

Table 2. Properties of the real model, the estimation model, and the simulation model. Errors in the simulation model and the estimation model are given in comparison to the real model. $s_{1_0},p_{p_0}$, and $p_{1_0}$ represent the initial actuator position, the initial pump pressure, and the initial pressure on the piston side as the system states. The system parameters ${\mathbf{C}}_{{\mathbf{v}}_a}$, ${\mathbf{C}}_{{\mathbf{v}}_b}$, ${\mathbf{C}}_{{\mathbf{v}}_c}$, ${\mathbf{C}}_{{\mathbf{v}}_d}$, $k_0$, and $k_p$ represent the semi-empiric flow rate coefficient at the a, b, c, and d ports of the directional control valve, the flow gain, and the pressure flow coefficients, respectively.

Errors	Symbol	Real Model	Estimation Model	Simulation Model
State	$s_{1_0}$	$\sqrt{3}\mathrm{m}$	$1.62\mathrm{m}$	$1.62\mathrm{m}$
State	$p_{p_0}$	7.6 MPa	5.6 MPa	5.6 MPa
State	$p_{1_0}$	1 MPa	2 MPa	2 MPa
Parameter	${\mathbf{C}}_{{\mathbf{v}}_a}$	Non-linear	Linear	Linear
Parameter	${\mathbf{C}}_{{\mathbf{v}}_b}$	Non-linear	Linear	Linear
Parameter	${\mathbf{C}}_{{\mathbf{v}}_c}$	Non-linear	Linear	Linear
Parameter	${\mathbf{C}}_{{\mathbf{v}}_d}$	Non-linear	Linear	Linear
Parameter	$k_0$	0.5	0.4	0.4
Parameter	$k_p$	1600 MPa	1500 MPa	1500 MPa

, the properties of the real model, the estimation model, and the simulation model are provided. Note that the simulation model is used in this study to demonstrate the differences between the simulated world and the real world.

As in practise, the minimum and maximum points on the characteristic curves of a directional control valve can be determined from the manufacturer’s catalogues. Using this limited information, the characteristic curves are defined linearly at all ports of the directional control valve in the cases of the estimation model and the simulation model. The linear characteristic curves are implemented by using the minimum and maximum values of the semi-empiric flow rate coefficients ${\mathbf{C}}_{{\mathbf{v}}_a}$, ${\mathbf{C}}_{{\mathbf{v}}_b}$, ${\mathbf{C}}_{{\mathbf{v}}_c}$, and ${\mathbf{C}}_{{\mathbf{v}}_d}$ at the valve closing and the valve opening positions, respectively. The linear characteristic curves of the directional control valve affect the dynamics of the estimation model throughout the simulation runtime. In the case of the real model, the characteristic curves of the directional control valve are unclear and can be non-linear. With Equation (21), the non-linear characteristic curves of the directional control valve are implemented using ${\mathbf{C}}_{{\mathbf{v}}_a}$, ${\mathbf{C}}_{{\mathbf{v}}_b}$, ${\mathbf{C}}_{{\mathbf{v}}_c}$, and ${\mathbf{C}}_{{\mathbf{v}}_d}$ in the hydraulic circuit of the real model. Similarly, the initial actuator positions $s_{1_0}$ of the real model and the estimation model are different. Note that the initial relative joint coordinates of the bodies in the system can be found from $s_{1_0}$ and ${\dot{s}}_{1_0}$ by using geometrical relationships. To avoid instabilities in the integration process, the simulations are started in the static equilibrium position, the details of the mechanism of which can be found in [34].

3.1.2. Sensor Measurements

In this study, the measurable observations $\mathbf{o}={\left[\begin{array}{ccc}s& p_p& p_1\end{array}\right]}^T$ are taken from the real model. In the real model, the actuator position sensor measures the actuator position s [76]. Gauge pressure sensors are used for the pressure measurements $p_p$ and $p_1$ [34]. These pressure sensors measure the pressure with respect to the atmospheric pressure. The pressure sensors are installed on their respective volumes, as also shown in Figure 2. The numerical values of the standard deviation, as mentioned in Equation (29), are taken as ${\sigma }_s^{\prime }=1.12\times {10}^{-3}\mathrm{m}$ and ${\sigma }_p^{\prime }=1.5\times {10}^5$ Pa for the actuator and pressure sensors, respectively.

3.2. Parameter Estimation Algorithm

In the parameter estimation algorithm, the augmented state vector ${\widehat{\mathbf{x}}}^{\prime }$ is defined as

(37)

where s is the actuator position, $\dot{s}$ is the actuator velocity, $p_p$, $p_1$, and $p_2$ are the pressures, $k_p$ is the pressure flow coefficient, $k_0$ is the flow gain, and ${\mathbf{C}}_{\mathbf{v}}=\left[\begin{array}{cccc}{\mathbf{C}}_{{\mathbf{v}}_a}& {\mathbf{C}}_{{\mathbf{v}}_b}& {\mathbf{C}}_{{\mathbf{v}}_c}& {\mathbf{C}}_{{\mathbf{v}}_d}\end{array}\right]$ are the semi-empiric flow rate coefficients at the corresponding ports of the directional control valve. In Equation (37), $\widehat{\mathbf{x}}$ and $\widehat{\mathbf{y}}$ present the states and the parameters of the hydraulically driven four-bar mechanism, respectively. Equations (20)–(28) are implemented to estimate the augmented state vector ${\widehat{\mathbf{x}}}^{\prime }$ and the characteristic curves of the directional control valve. In this application, the third-order B-spline interpolation method is combined with the ADEKF. For the case example, three, four, five, and six control points are used in the parameter estimation algorithm to compute Equations (20) and (24). As mentioned earlier, the first, third, and fourth state variables are measured. Therefore, the sensor measurement function $\mathbf{h}\left({\widehat{\mathbf{x}}}_k^{\prime -}\right)$ and its Jacobian ${\mathbf{h}}_{{\mathbf{x}}^{\prime }}$ can be written as

$$\left.\begin{array}{cc}\hfill \mathbf{h}\left({\widehat{\mathbf{x}}}_k^{\prime -}\right)& ={\left[\begin{array}{ccc}{\widehat{\mathbf{x}}}_{k,1}^{\prime -}& {\widehat{\mathbf{x}}}_{k,3}^{\prime -}& {\widehat{\mathbf{x}}}_{k,4}^{\prime -}\end{array}\right]}^T\hfill \\ \hfill {\mathbf{h}}_{{\mathbf{x}}^{\prime }}& =\left[\begin{array}{ccccccccccc}1& 0& 0& 0& 0& 0& 0& {\displaystyle \frac{\partial \mathbf{h}}{\partial {\mathbf{C}}_{{\mathbf{v}}_a}}}& {\displaystyle \frac{\partial \mathbf{h}}{\partial {\mathbf{C}}_{{\mathbf{v}}_b}}}& {\displaystyle \frac{\partial \mathbf{h}}{\partial {\mathbf{C}}_{{\mathbf{v}}_c}}}& {\displaystyle \frac{\partial \mathbf{h}}{\partial {\mathbf{C}}_{{\mathbf{v}}_d}}}\\ 0& 0& 1& 0& 0& 0& 0& {\displaystyle \frac{\partial \mathbf{h}}{\partial {\mathbf{C}}_{{\mathbf{v}}_a}}}& {\displaystyle \frac{\partial \mathbf{h}}{\partial {\mathbf{C}}_{{\mathbf{v}}_b}}}& {\displaystyle \frac{\partial \mathbf{h}}{\partial {\mathbf{C}}_{{\mathbf{v}}_c}}}& {\displaystyle \frac{\partial \mathbf{h}}{\partial {\mathbf{C}}_{{\mathbf{v}}_d}}}\\ 0& 0& 0& 1& 0& 0& 0& {\displaystyle \frac{\partial \mathbf{h}}{\partial {\mathbf{C}}_{{\mathbf{v}}_a}}}& {\displaystyle \frac{\partial \mathbf{h}}{\partial {\mathbf{C}}_{{\mathbf{v}}_b}}}& {\displaystyle \frac{\partial \mathbf{h}}{\partial {\mathbf{C}}_{{\mathbf{v}}_c}}}& {\displaystyle \frac{\partial \mathbf{h}}{\partial {\mathbf{C}}_{{\mathbf{v}}_d}}}\end{array}\right]\hfill \end{array}\right\},$$

(38)

where ${\mathbf{C}}_{\mathbf{v}}=\left[\begin{array}{cccc}{\mathbf{C}}_{{\mathbf{v}}_a}& {\mathbf{C}}_{{\mathbf{v}}_b}& {\mathbf{C}}_{{\mathbf{v}}_c}& {\mathbf{C}}_{{\mathbf{v}}_d}\end{array}\right]$. $\mathbf{h}\left({\widehat{\mathbf{x}}}_k^{\prime -}\right)$ and ${\mathbf{h}}_{{\mathbf{x}}^{\prime }}$ are used in Equations (26) and (27), respectively, in the parameter estimation algorithm.

4. Results and Discussion

In this section, the results of the simulation of the estimation model of the parameter estimation algorithm are presented. The results of the estimation model are compared to those of the real model and the simulation model. The initial covariance ${\mathbf{P}}_0$ used in the augmented state estimator includes ${\sigma }_s^2=1\times {10}^{-4}{\mathrm{m}}^2$ for the actuator position, ${\sigma }_{\dot{s}}^2=1\times {10}^{-4}{\mathrm{m}}^2/{\mathrm{s}}^2$ for the actuator velocity, and three pressure terms of ${\sigma }_p^2=22.50\times {10}^7{\mathrm{Pa}}^2$ in the diagonal. For the hydraulic parameters, the initial covariance values ${\sigma }_{k_p}^2=1\times {10}^{14}{\mathrm{Pa}}^2$, ${\sigma }_{k_0}^2=1\times {10}^2$ and ${\sigma }_{{\mathbf{C}}_{\mathbf{v}}}^2=9\times {10}^2\frac{m^6}{s^{2}Pa}$ are used in the diagonal. The numerical values of the plant noise ${\sigma }_{\ddot{\mathbf{x}}}^2=0.8{\mathrm{m}}^2/{\mathrm{s}}^4$ and ${\sigma }_p^2=259.81\times {10}^7{\mathrm{Pa}}^2$ for Equation (31) are obtained through trial and error. All models are run with a time step of $1$m$\mathrm{s}$ and provide sensor data to the parameter estimation algorithm at 1000 Hz.

4.1. Estimating the Characteristic Curve of the Valve

In the real model, as only the minimum point ${\mathrm{c}}_{\min }$ and the maximum point ${\mathrm{c}}_{\max }$ on the characteristic curves are known at the a, b, c, and d ports of the directional control valve, the characteristic curves are generally unclear in the working cycles of the real model. The characteristic curves may vary from one valve to another and can be highly non-linear in the working cycle. In Figure 3, Spline 1 and Spline 2, which are in the cyan colour, are used to demonstrate the non-linear behaviour of the directional control valve in the real model.

The proposed parameter estimation algorithm can be used to estimate the characteristic curves of the real model with this limited information. To this end, the semi-empiric flow rate coefficients ${\mathbf{C}}_{{\mathbf{v}}_a}$, ${\mathbf{C}}_{{\mathbf{v}}_b}$, ${\mathbf{C}}_{{\mathbf{v}}_c}$, and ${\mathbf{C}}_{{\mathbf{v}}_d}$ are defined with the data control points between ${\mathrm{c}}_{\min }$ and ${\mathrm{c}}_{\max }$ in Equation (37). Equation (37) is further used in terms of the control point vector $\mathbf{N}$ in Equations (20)–(28) to estimate the characteristic curves.

For instance, in the case of Figure 4, the control point vector ${\mathbf{N}}_a$ at port a of the directional control valve can be defined in terms of ${\mathrm{c}}_1$, ${\mathrm{c}}_2$, ${\mathrm{c}}_3$, and ${\mathrm{c}}_4$ as

$${\mathbf{N}}_a={\left[\begin{array}{cccccc}{\mathrm{c}}_{\min }& {\mathrm{c}}_1& {\mathrm{c}}_2& {\mathrm{c}}_3& {\mathrm{c}}_4& {\mathrm{c}}_{\max }\end{array}\right]}^T.$$

(39)

To estimate the characteristic curve, three, four, five, and six control points are used in the control point vector ${\mathbf{N}}_a$. As an example, these data control points for Spline 1 in each case are presented in

Table 3. Root mean square error in the estimation of the characteristic curve. The third and fourth columns represent the root mean square errors in Spline 1 and Spline 2, respectively.

Control Points	Control Point Vector ${\mathbf{N}}_{\mathit{a}}$	RMSE	RMSE
Three points	$\left[\begin{array}{ccc}0& 5& 10\\ 0& 47.5& 95\end{array}\right]$	$0.04\%$	$0.03\%$
Four points	$\left[\begin{array}{cccc}0& 3.3& 6.6& 10\\ 0& 31.6& 63.3& 95\end{array}\right]$	$0.05\%$	$0.01\%$
Five points	$\left[\begin{array}{ccccc}0& 2.5& 5& 7.5& 10\\ 0& 23.7& 47.5& 71.2& 95\end{array}\right]$	$0.06\%$	$0.07\%$
Six points	$\left[\begin{array}{cccccc}0& 2& 4& 6& 8& 10\\ 0& 19& 38& 57& 76& 95\end{array}\right]$	$0.07\%$	$0.08\%$

.

The abscissa of vector ${\mathbf{N}}_a$ represents the spool position U, whereas the ordinate of vector ${\mathbf{N}}_a$ indicates the semi-empiric flow rate coefficients ${\mathbf{C}}_{{\mathbf{v}}_a}$ at port a of the directional control valve. The results of the second-order B-spline are described in Appendix A. The results of the third-order B-spline demonstrate the characteristic curve of the directional control valve relatively better in a working cycle, as shown in Figure 3. As can be seen, Spline 1 and Spline 2 of the third-order B-spline are drawn in each data control point estimation case. The dashed red-coloured line indicates the characteristic curve of the simulation model. The dashed black-coloured line demonstrates the estimation model. In Figure 3a, three points, ${\mathrm{c}}_{\min }$, ${\mathrm{c}}_1$, and ${\mathrm{c}}_{\max }$, are used to estimate Spline1 and Spline2 of the real model. The characteristic curve of the estimation model precisely follows the real model in the case of three points. Further, the percentages of the root mean square error (RMSE) are described in the Table 3 for Spline 1 and Spline 2 to verify the observations.

Figure 3b shows the estimation of the characteristic curve when using four points ${\mathrm{c}}_{\min }$, ${\mathrm{c}}_1$, ${\mathrm{c}}_2$, and ${\mathrm{c}}_{\max }$. As can be seen in Figure 3b, the semi-empiric flow rate coefficient ${\mathbf{C}}_{{\mathbf{v}}_a}$ for Spline 2 changes with small increments until $52\%$ opening of the spool as compared to Spline1 in the real model. After this point, the parameter ${\mathbf{C}}_{{\mathbf{v}}_a}$ increases sharply towards the maximum point ${\mathrm{c}}_{\max }$. The difference of the estimated curve from the real model’s curve is indistinguishable. The RMSEs of these curves are given in Table 3. The relatively complicated non-linear behaviours of the directional control valve can be estimated by using five control points and six control points. This can be seen in Figure 3c,d. By using the estimated characteristic curves, the working conditions of the directional control valve can be predicted.

4.2. Convergence of the Vector Data Control Points

The convergence rate of the data control points in the parameter vector ${\mathbf{C}}_{{\mathbf{v}}_a}$ is further explained in Figure 5 to describe the estimation process. These plots demonstrate the convergence rate of data control points in the case of Spline 2, as presented in Figure 3. For instance (see Figure 5b), ${\mathrm{c}}_1$ and ${\mathrm{c}}_2$ converge towards the corresponding point on the curve of the real model at $0.22\mathrm{s}$. However, during the estimation process, ${\mathrm{c}}_1$ briefly becomes negative, and shortly thereafter converges smoothly to the real model.

The curves of Spline 2 change into an S-shape during the working cycle, as shown in Figure 3c,d. In these cases, ${\mathrm{c}}_2$, ${\mathrm{c}}_3$, and ${\mathrm{c}}_4$ converge at different simulation times according to the corresponding order in the vector ${\mathbf{C}}_{{\mathbf{v}}_a}$. Through the ADEKF algorithm, unknown curves start converging within a range of $0<t\le 0.3\mathrm{s}$ when using the three-, four-, five-, and six-point estimation techniques.

4.3. Accuracy Requirements of State Estimations

The successful application of the parameter estimation algorithm requires the accurate estimation of the system states $\mathbf{x}$. To demonstrate this requirement, the errors in the estimated actuator position s, estimated actuator velocity $\dot{s}$, estimated pump pressure $p_p$, estimated piston side pressure $p_1$, estimated piston-rod side pressure $p_2$, and estimated parameter ${\mathbf{C}}_{{\mathbf{v}}_a}$ in the case of Spline 2 (described in Figure 3d) are shown in Figure 6. The errors in the estimated parameters $k_p$ and $k_0$ are presented in Appendix B. The average of the parameter vector ${\mathbf{C}}_{{\mathbf{v}}_a}$ at each time step is considered in Figure 6.

The errors are computed from $\pm 1.96\sigma$. Here, $\sigma$ is the standard deviation calculated from the covariance matrix ${\mathbf{P}}_k^{+}$ at each time step. These plots demonstrate the requirement of an accurate estimation of the system’s states to estimate the system’s parameters. As can be seen in Figure 6, the $95\%$ confidence interval (CI) is used by the system states in the $5\mathrm{s}$ simulation period. The errors in s, $\dot{s}$, $p_p$, $p_1$, and $\dot{s}$ fluctuate in the confidence interval. As indicated earlier, s, $\dot{s}$, and $p_p$ are measured in this example. The errors in the parameters ${\mathbf{C}}_{{\mathbf{v}}_a}$, $k_p$, and $k_0$ are also in the CI, as can be seen in the corresponding plots. The key to the parameter estimation is that the estimated system states should be in the $95\%$ CI during the working cycle.

5. Conclusions

This work proposes the estimation of the parameters of a system by combining parameter estimation theories and curve-fitting methods. The ADEKF algorithm is introduced in the framework of a B-spline curve-fitting method. Using the proposed algorithm, the parameters can be defined as a vector containing a set of data control points. This algorithm is applied on a hydraulically driven four-bar mechanism to estimate the characteristic curves of a directional control valve. The double-step semi-recursive formulation and lumped fluid theories are used to model the four-bar mechanism and the hydraulic system, respectively. The measurements taken from the real system include the actuator position, pump pressure, and piston side pressure. The semi-empiric flow rate coefficient vector ${\mathbf{C}}_{{\mathbf{v}}_a}$ is defined with three to six data control points in order to define the characteristic curve of the directional control valve.

The unknown non-linear nature of the characteristic curves of the directional control valve are precisely estimated. The maximum RMSE observed in the estimation of the characteristic curves is 0.08%. This implies that the characteristic curves are accurately estimated. The data control points in the parameter vector ${\mathbf{C}}_{{\mathbf{v}}_a}$ converge in the range of $0<t\le 0.3\mathrm{s}$ in these estimation cases. To account for the system’s response, the estimation of the system’s state vector variables should be located in the 95% confidence interval. By using the estimated characteristic curves, important information about the discharge coefficient, pressure losses, and flow characteristics of the directional control valve can be interpreted. With this valuable information, manufacturers and users can monitor the condition of a system and make decisions about the repair and maintenance of hydraulically driven systems.

Applying the parameter estimation algorithm in the real world by using a multibody-based estimation model can enable the estimation of important parameters. This can be challenging, as the estimation model might not be as accurate as the real world necessitates. However, despite implementation challenges, the application of this parameter estimation algorithm will provide an interesting area for manufacturers and researchers. Manufacturers can use these parameters in condition monitoring, repair and maintenance, and the anticipation of product life cycles. With a product’s application history, important design changes can be introduced in future designs of the product. This will ultimately lead to more efficient MBS-based digital-twin applications through the use of real-time simulations and more sustainable future products.