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 [21, 22]. Parameter estimation through
3.2 Estimation algorithm 53
Errors in the model
Sensors
h model
Simulation (multibody and
hydraulics)
ok
Inputs
ADEKF ˆ
x0−k ,P−k
xˆ0+k ,P+k ˆ
x0−k−1, P−k−1
states and Estimated parameters
Figure 3.1. Parameter estimation methodology.
the ADEKF comprises prediction and update stages. At the prediction stage, in case of the coupled multibody and hydraulic systems, the augmented state vector can be described asx0=h(zi)T ( ˙zi)T pT yTiT. At this step, ˆx0−k is calculated in the time integration of a dynamic model described as [68]
ˆ
x0−k =f(ˆx0+k−1, Uk). (3.1) To account for unknown parameters, the proposed parameter estimation algorithm employs the curve-fitting method. Through this method, the parameters are introduced in Eq. (3.1) as the knot vectoru. A B-spline curveC(u) is constructed for non-uniform open splines [21, 22] at the current time step to represent parameters as,
C(u) =
n
X
i=0
Bi,d(u)Ni, (3.2)
where n is the number of control points,i is theith number of control points,dis the degree, Bi,d(u) are thedth order of B-spline basis functions, and Ni is the control point vector. The control point vector can be expressed in terms of the system parameters y. For instance, in case of the characteristic curve, the control point vector can be written in terms of the spool position and semi-empiric flow
54 3 Parameter estimation method
rate coefficient as N=
"
Umin U1 ... Un Umax
Cvmin Cv1 ... Cvn Cvmax
#
. Here, Umin, U1, and Un represent spool positions, andCvmin, Cv1, andCvmax 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 [21,22] as,
Bi,0(u) =
1 if ui≤u<ui+1
0, otherwise , (3.3)
Bi,j(u) = u−ui
ui+j−uiBi,j−1(u) + ui+j+1−u
ui+j+1−ui+1Bi+1,j−1(u), (3.4) whereui is theith element of the knot vector for non-uniform open splines. Next, the numerical values of parameters, which are scalar, should be evaluated by using Eq. (3.2) at time step k to be incorporated in Eq. (3.1). The calculation of Eq. (3.1) at the desired input signal is straightforward. However, the numerical computation of the Jacobian matrix fx0
k−1 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 [43, 71] for very small increments. For instance, in the case of a multi-variable function, the Jacobian column of Eq. (3.1) with respect to therth term of the augmented state vector ˆ
x0k−1 can be written in the partial derivative form as,
∂f(ˆxk−1,10 , ...,xˆ0k−1,L)
∂xˆ0k−1,r = Im(f(ˆx0+k−1,1, ...,ˆx0k−1,r+iδ, ...,ˆx0+k−1,L))
δ +O(δ2), (3.5)
wherer∈[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 rth term of the state vector ˆx0k−1,r+iδ is expanded by using the Taylor series [43, 71].
The evaluation of Im(f(ˆx0+k−1,1,xˆ0+k−1,2, ...,xˆ0k−1,r + iδ, ...,xˆ0+k−1,L)) requires the calculation of C(uk−1) as complex numbers to include the parameter vector yk−1. However, the knot vector cannot contain any complex numbers [21, 22].
Therefore, a criterion|u−uik−1|< Ξ is introduced, whereΞ can be a small real number, such as 0.1, which implies that the knot-point distance between uik−1 and the complex input argument ushould be less thanΞ. Using this criterion enables the knot points to be evaluated with the curve-fitting method through
3.2 Estimation algorithm 55
complex input. With the Jacobian of the dynamic systemfx0
k−1, the covariance matrix P−k is approximated as [68],
P−k =fx0
k−1P+k−1fxT0
k−1+Qk, (3.6)
whereQk 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 ˆxk0−. The innovation∆k is calculated as [68],
∆k =ok−h(ˆx0−k ), (3.7)
whereok are the sensor measurements at the kth time step, and h(x0−k ) is the sensor measurement function. With the Jacobian of the sensor measurement functionhx0
k, the innovation in the covariance matrixSk and the Kalman gain Kk can be calculated as [68],
Sk =hx0
kP−khTx0 k+Rk Kk =P−khxT0
kS−1k
, (3.8)
where Rk is the covariance matrix of the measurement noise. Finally, the augmented state vector ˆx0+k and covariance matrixP+k are updated at the time stepk, which will be used for the next time step as [68],
ˆ
x0+k = ˆx0−k +Kk∆k
P+k = (IL−Kkhx0
k)P−k
, (3.9)
whereIL is the identity matrix of dimensionL.
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 mea-surement noises. Furthermore, it was established in [64,63] 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.
56 3 Parameter estimation method
For instance, when using position and pressure sensors, the covariance matrix of the measurement noise,R, would then take the form [37, 63, 64],
R=
whereσs0 andσ0pare the standard deviations of measurement noises at the position and pressure levels, respectively. In Eq. (3.10), n is the number of actuator sensors andnp is the number of pressure sensors. In, Inp, 0n×np, and 0np×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 [64, 63] and can be written as,
where∆tis the size of the integration time step andnf is the number of degrees of freedom of the system. It should be noted that Eq. (3.11) 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 [37], the variance of hydraulic pressures,σpD, and the variance of hydraulic parameters,σhpD, can be directly incorporated as the diagonal elements in Eq. (3.11). Accordingly, the structure of the covariance matrix of the plant noise, Q, in the parameter estimation can be written as
Q=
In this study, the integration errors are assumed to be negligible in comparison to the acceleration, pressure, and parameter errors.
Chapter 4
Summary of findings
The main motivation of this dissertation is to develop methods to integrate multibody-dynamics-driven digital solutions into the product lifecycle, and in-vestigate the technical aspects associated to the user experiences in real-time simulation applications. This chapter describes important findings of the used publications to achieve the objective of this dissertation.
4.1 Real-time simulation monolithic approaches
Using monolithic approaches, the hydraulic systems can be combined with a multibody approach in the framework of real-time simulation. The internal dynamics of the hydraulics is defined with the resultant hydraulic force in terms of actuator position and actuator velocity. The modeling of hydraulic systems results in numerical stiffness. This problem can be reduced by proper selection of a multibody approach. Naya et al. [53] combined the index-3 augmented Lagrangian formulation with the hydraulics in the monolithic studies using the global coordinate system. In turn, Rahikainen et al. implemented monolithic approaches based on the relative coordinates for the coupled simulation of multibody and hydraulic dynamics [57,58]. Publication I introduces a monolithic scheme for the coupled simulation of the double-step semi-recursive formulation and hydraulic systems.
Publication I
In this study, the index-3 augmented Lagrangian and double-step based semi-recursive formulations with hydraulic systems are compared. As described in chapter 2, the double-step semi-recursive formulation employs a coordinate
57
58 4 Summary of findings
partitioning approach [60] whereas, the index-3 augmented Lagrangian approach uses a full set of coordinates [6,15] to enforce constraint equations. The hydraulics systems are modelled using lumped fluid theory as described in section 2.2, and coupled with the multibody equations of motion according to section 2.3 in a monolithic study.
InPublication I, the hydraulically actuated quick-return and four-bar mechanisms are used as numerical examples. The modeling details of hydraulic systems, initial conditions and the physical parameters of the system can be found in section 5 of Publication I. The results of the double-step and the index-3 augmented Lagrangian semi-recursive formulations are compared based on the constraint violation, numerical efficiency, energy balance and work cycle.
As can be seen in Figure 4.1, the quick-return mechanism includes two closed loops.
However, in Publication I, the results of only one loop are presented. In both mechanisms, the constraints are fulfilled with good accuracy at position, velocity and acceleration levels with the multibody approaches. See Figure 11 and Figure 12 inPublication I for further details. The double-step formulation offers advantages by fulfilling the constraints to the machine precision level and maintaining good error control. However, the index-3 augmented Lagrangian semi-recursive formulation can provide accurate solutions to the singular configurations [27] and redundant constraints.
Joints at points O, H, I, K, andL: Revolute joint Joints at points J and M: Translational joint (cut-joint)
Q
M
p3, Be3, V3
Z
Figure 4.1. A quick-return mechanism actuated by a double-acting hydraulic cylinder used inPublication I.