Determining the State of a Nonlinear Flexible Multibody System using an Unscented Kalman Filter

This paper describes an estimator incorporating the Unscented Kalman Filter (UKF) technique and multibody system dynamics, to determine the state of the flexible multibody applications. The dynamic equation of the flexible mechanism is formed using a set of non-linear equations as functions of reference and modal coordinates. Since both reference and modal coordinates have no physical meaning, their information is not able to be obtained directly from sensors. Thus, a novel technique is proposed in this work that can successfully translate physical measurements collected by sensors into non-physical modal coordinates. To validate the performance of the proposed modeling technique to apply a UKF to determine the state of a nonlinear flexible multibody system, simulation were carried out for a four-bar mechanism case study to compare the simulation data and UKF data.

The global position of an arbitrary particle P within a 192 deformable body i in the FFRF, as shown in Fig. 1 can be 193 written as: where R is the position vector of the floating frame of the displacement vectorū P f is given by: where Φ t is the translational modal transformation matrix where ϕ is the orientation of the floating frame of reference.

211
Equations of motion can be written using independent 212 coordinate systems, which makes it possible to convert the 213 system to state space forms [34]. According to the structure 214 of the dynamic system and the number of degrees of freedom, 215 the coordinate of a mechanical system, q, can be partitioned 216 into dependent and independent coordinates, q d and q i , re-217 spectively, as: where the number of dependent coordinates is equal to the where C d q and C i q are the partitioned part of the constraint 229 Jacobian matrix related to the dependent and independent 230 coordinates, respectively.

231
Subsequently, the virtual change of the mechanical system 232 coordinates can be written as: where matrix B i is given by: where I is the identity matrix with the size of n c × n c .

235
Thereafter, the second derivative ofq can be written as: where matrix D is expressed as follows, where 0 is a null matrix with the size of n c × 1.

238
The equations of motion can be formulated according to where M is the mass matrix, C q is the Jacobian matrix of 245 constraint equations, λ is the vector of Lagrange multipliers,

246
Q v is the quadratic velocity vector and Q e is the external 247 force vector, which comprises of external and elastic forces.

248
Substituting Eqs. (7) and (9) into Eq. (11) results in: where the virtual velocities δq i includes the independent velocities δq i are independent the following can be written: In in whichẋ is the time derivative of x and f (x, t) represents 275 the dynamic system function. estimate the data (x k , and P k ). Furthermore, the estimated 294 data (x k , and the error covariance matrix P k ) feeds back into 295 the Kalman filter for the next step.

296
Note that the measurements are taken from real system. In qï where w is the process noise that occurs during the simulation of the system. Moreover, measurement data must be defined through an equation to represent the linear relationship between the sensors and system states [34]. Therefore, the output function of the system y can be expressed as: where h is the measurement sensitivity matrix, a linear dif-  The next step is to discretize the continuous-time system given in Eq. (16). This is done using the Runge kutta fourthorder method. Therefore, for the Eq. (16) and Eq. (17), the discontinuous system can be expressed as: where k is the time step. F and H symbolize the discontinu- points. Based on the assumptions ofx 0|0 and P 0|0 the sigma 329 points are calculated for k = 1 as follows: where L is the number of states and j = 0, ..., 2L, 331 P k−1|k−1 is the matrix square root of the Cholesky decom-332 position using its lower triangular matrix, Γ = α 2 (L + κ).

333
Scaling factors κ control the weighting of the sigma points.
Thereafter, the mean and covariance should be calculated.

342
The predicted covariance and mean are calculated based on 343 the results of the sigma points as follows, where w , β is a secondary scaling factor.

346
The value of β is considered to be 2 for Gaussian noise 347 distributions.

348
Next, predicted measurement points are computed with the predicted state points γ (j) k|k−1 . At this stage, determining the prediction state variables is accomplished by propagating sigma points in the linear measurement sensitivity matrix, The predicted measurement i.e. the mean of the predicted measurement points, (ŷ), is calculated next. Specifically, these are the predicted state variables with the scaling factor, The innovation covariance and the cross covariance for the predicted measurement can then be computed as: where R k is the covariance for the measurement noise. To rectify the states and the covariance matrix, the Kalman gain matrix should be calculated [2], This is then used to calculate the a posteriori estimated state and covariance matrix:  Fig. 3 are related as follows, whereĨ 2 is a skew symmetric matrix as    As the figure shows, the modal matrix is partitioned into 399 the master and slave matrices, Φ ma t , and Φ sl t , respectively, 400 Therefore, the nodal coordinates can be written in two parts.

401
Subsequently, following equations can be written: Thereafter the modal coordinates can be written as follows, The modal matrix, Φ ma t can be partitioned according to the consists of four elements. Table 1 gives the specifications for 461 each beam of the four-bar mechanism.

462
The boundary condition between each body is the revolute    The node next to and after the middle node was also stud-517 ied with respect to different initial conditions. Fig. 9 shows