Online Dynamic Load Identification Based on Extended Kalman Filter for Structures with Varying Parameters

Abstract

Dynamic load identification is an inverse problem concerned with finding the load applied on a structure when the dynamic characteristics and the response of the structure are known. In engineering applications, some of the structure parameters such as the mass or the stiffness may be unknown and/or may change in time. In this paper, an online dynamic load identification algorithm based on an extended Kalman filter is proposed. The algorithm not only identifies the load by measuring the structural response but also identifies the unknown structure parameters and tracks their changes. We discuss the proposed algorithm for the cases when the unknown parameters are the stiffness or the mass coefficients. Furthermore, for a system with many degrees of freedom and to achieve online computations, we implement the model reduction theory. Thus, we reduce the number of degrees of freedom in the resulting symmetric system before applying the proposed extended Kalman filter algorithm. The algorithm is used to recover the dynamic loads in three numerical examples. It is also used to identify the dynamic load in a lab experiment for a structure with varying parameters. The simulations and the experimental results show that the proposed algorithm is effective and can simultaneously identify the parameters and any changes in them as well as the applied dynamic load.

1. Introduction

In many applications related to structural safety and reliability, it is necessary to obtain the dynamic loads active on a structure. Furthermore, it is possible that the dynamic parameters of such a structure vary with time. In such cases, especially for gradually varying parameters, the dynamic characteristics of a structure cannot be known in advance. Therefore, in order to identify the load acting on such structures, it is important to identify the structural parameters before or simultaneously as we identify the load by solving an inverse problem starting from the structural responses. Such information will not only be useful for monitoring the loads but also to monitor the characteristics of a structure as it changes in real time.

The methods of parameters’ identification can be divided into time-domain and frequency-domain methods [1,2,3,4]. In the frequency-domain, the methods are in general concerned with identifying the modal parameters of a structure, and then using the modal coordinate transformation relationship to obtain the physical parameters [5,6,7]. As for the time-domain, several scholars have proposed a variety of identification methods for the modal parameters as well as the physical parameters. The proposed methods include: Random Decrement Technique (RDT) [8], Ibrahim Time Domain Method (ITD) [9], Stochastic Subspace Identification (SSI) [10], Eigensystem Realization Algorithm (ERA) [11] and Natural Excitation Technique (NExT) [12,13]. However, using these methods requires the excitation to be white noise, or the response to be a free decay process. There are also extended Kalman filter (EKF) [14], least Square Estimation (LSE) [15], H$\infty$filter [16] and Monte Carlo filter [17] to identify the physical parameters for a structure in the time-domain.

Similarly, the dynamic load identification can also be classified into frequency- and time-domain methods. The frequency-domain methods are relatively more mature, but they have certain limitations. The identification accuracy strongly depends on the time length of the response data acquisition. Their performance is negatively impacted if the response data are only available for a short time period. The time-domain methods are relatively recent compared to the frequency-domain. In general, time-domain methods express the system characteristics as a Duhamel integral, and then obtain the load acting on a structure using a convolution response. Eills [18] is the first to identify the aerodynamic load acting on a structure using the frequency-domain. Bartlett [19] identified the load on the hub by measuring the acceleration response of a helicopter. Inoue [20] used deconvolution method to identify the load in the frequency-domain. Milana [21] proposed a method of load identification by coherent analysis of a structural response. Liu et al. [22] used an interpolation function to fit the load to be identified. Zhang [23] deduced the method of identifying two dimensional distributed dynamic loads in generalized orthogonal domains. By introducing the Tikhonov regularization algorithm, Jiang [24,25] improved the accuracy of identifying distributed random dynamic loads in the generalized orthogonal polynomial domain. Hory [26] first used the method of time discretization to identify aircraft flight loads. Wu et al. [27] identified the load in time domain based on Green’s function. Gunawan [28] improved the recognition effect of the impact load with the help of a L-curve and the truncated singular value method. Song [29] proposed a dynamic load identification method based on the Kalman filter, which took the strain signal as the observation signal and identified the unknown load on a cantilever beam through recursive iteration. The above methods can only identify the excitation when the dynamic characteristics of the structure are known.

Starting in the past decade, a lot of efforts were invested in solving the inverse problem of identifying the loads when the structural parameters are unknown while the dynamic model is given. Based on the classical Kalman filter (KF) [30], Yang [31] proposed an adaptive Kalman filtering algorithm. T. Lenkutis [32] modified the sinusoidal approximation using a windowing function, which can remove sharp jumps in a generated road profile. Y. Lei et al. [33] used an extended Kalman filter to identify the load and the state vectors through multiple stages. The algorithm can identify the dynamic load when the structural parameters are unknown and track changes in the structural parameters. Naets et al. [34] reconstituted the structural state vector parameters and excitation into an augmented vector, and then used the extended Kalman filter to identify the load and the unknown parameters. Du et al. [35] studied the identification of the blade tip vortex and the motion characteristics of the vortex. However, the “Drift phenomenon” is observed in the displacement and in the load identification if the measured acceleration is contaminated with noise. To solve this problem, Huang et al. [36] proposed a data fusion method, which used partial data fusion of displacement and acceleration responses in the observation equation. Lei et al. [37] proposed a dynamic load identification algorithm based on the recursive least square method. Zhang and He [38] introduced the projection matrix to identify jointly structural parameters and unknown loads. Yu [39] combined the theory of the Hilbert transform and Kalman filter to identify the parameters and the dynamic load of a shear frame structure.

In this paper, we are concerned with the dynamic load identification where structural parameters are unknown and change gradually in time. An example of such a system is a launched rocket where the mass is reduced in time due to the fuel consumption as well as moving through stages of a multistage rocket. Standard load identification methods cannot identify the load acting on such systems. We propose a load identification algorithm based on an extended Kalman filter. We validate the algorithm using numerical and experimental studies where we also consider data contaminated with noise. The algorithm is used to identify properties of structures with gradually varying stiffness or with gradually varying mass. To deal with continuous structures, they are first discretized using the finite element method. The discretization produces a symmetric linear system of equations. If the system includes many degrees of freedom, we implement the model reduction theory to reduce the number of degrees of freedom before applying the proposed extended Kalman filter algorithm. The rest of this paper is organized as follows: the load identification method is presented in Section 2; in Section 3, we discuss implementing the model reduction theory to reduce the number of degrees of freedom; in Section 4, three numerical examples are presented to study the efficiency and the accuracy of the proposed algorithm using different structural parameters and different noise conditions, while in Section 5, the proposed algorithm is used to identify the dynamic load of a simply supported beam with gradually varying mass characteristic in an experimental setting; we finish in Section 6 with some concluding remarks.

2. Inverse Algorithm

The general vibration differential equation of a n-DOF system is written as

$$M\ddot{p}\left(t\right)+C\dot{p}\left(t\right)+Kp\left(t\right)=B_{u}f\left(t\right).$$

(1)

Here, the parameters $M, C$ and $K$ represent the mass, the damping and the stiffness symmetric matrices of n × n for a given structure, while the variables $\ddot{p}(t)$, $\dot{p}\left(t\right)$ and $p\left(t\right)$ are acceleration, velocity and displacement vectors of n × 1, respectively.

Starting form Equation (1) and employing $x\left(t\right)=\left[\begin{array}{c}p\left(t\right)\\ \dot{p}\left(t\right)\end{array}\right]$ to describe the state vector of the structure, a linear state space of the vibration differential equation can be obtained with

$$\dot{x}\left(t\right)=\left[\begin{array}{cc}\mathbf{0}& I\\ -M^{-1}K& -M^{-1}C\end{array}\right]x\left(t\right)+\left[\begin{array}{c}\mathbf{0}\\ M^{-1}{B}_u\end{array}\right]f\left(t\right) .$$

(2)

where the influence matrix $B_u$ is the n × s matrix which is composed of 0 and 1 to represent the location of the applied excitation so that the position of the excitation on the structure is allocated 1 and otherwise is allocated 0. The matrix $\mathbf{0}$ is a zero matrix. The variable $f\left(t\right)$ is the excitation vector of s × 1.

If some of the parameters in the mass or the stiffness or the damping matrices are unknown, we combine these unknown parameters with the state vector $x\left(t\right)$ to form the extended state vector $z(t)$ in Equation (3)

$$z\left(t\right)=\left[\begin{array}{c}x\left(t\right)\\ \alpha \end{array}\right]=\left[\begin{array}{c}p\left(t\right)\\ \begin{array}{c}\dot{p}\left(t\right)\\ \alpha \end{array}\end{array}\right].$$

(3)

Therefore, the state updating equation of the extended state vector can be rewritten as

$$\dot{z}\left(t\right)=\left[\begin{array}{c}\dot{p}\left(t\right)\\ -M^{-1}Kp\left(t\right)-M^{-1}C\dot{p}\left(t\right)+M^{-1}{B}_{u}f\left(t\right)\\ \mathbf{0}\end{array}\right]=f^c\left(z\left(t\right),f\left(t\right)\right).$$

(4)

Clearly, the state updating equation with the extended state vector is a nonlinear equation.

The form of measurement available in the measurement updating equation is related to the type of observation that can be gathered from the structure. Usually, the acceleration values measured at the degrees of freedom are the easiest to get. The measurement updated equation can be written as

$$y\left(t\right)=H_0{M}^{-1}\left(-Kp\left(t\right)-C\dot{p}\left(t\right)+B_{u}f\left(t\right)\right)=h\left(z\left(t\right),f\left(t\right)\right).$$

(5)

Thus, the state updating equation and the measurement updating equation are nonlinear equations with respect to two variables. Applying the Taylor expansion, the above two nonlinear equations are transformed into a sum of polynomials. We retain only the first-order polynomials which results into the following

$$f^c\left(z\left(t\right),f\left(t\right)\right)\approx f^c\left(z_{k-1|k-1},{\widehat{f}}_{k-1}\right)+{\nabla }_z{f}_{k-1}^c·\left(z\left(t\right)-z_{k-1|k-1}\right)+{\nabla }_f{f}_{k-1}^c·\left(f\left(t\right)-{\widehat{f}}_{k-1}\right),$$

(6)

$$h\left(z\left(t\right),f\left(t\right)\right)\approx h\left(z_{k|k-1},{\widehat{f}}_{k-1}\right)+{\nabla }_z{h}_k·\left(z\left(t\right)-z_{k|k-1}\right)+{\nabla }_f{h}_k·\left(f\left(t\right)-{\widehat{f}}_{k-1}\right).$$

(7)

where ${\nabla }_z{f}_{k-1}^c$ and ${\nabla }_z{h}_k$ are the first-order partial derivatives of the state updating equation and the measurement updating equation to the extended state vector $z\left(t\right)$, respectively. Then, ${\nabla }_f{f}_{k-1}^c$ and ${\nabla }_f{h}_k$ are the first-order partial derivatives of the state updating equation and the measurement updating equation to the load vector. The subscript $k$ or $k-1$ denotes the time step. The parameter $z_{k|k-1}$ indicates the prior estimation of the extended state vector at the $k$th time step starting from the previous time step, namely, $k-1$.

In order to identify the extended state vector and the unknown load, the extended Kalman filter can be divided into three steps: state updating, excitation identification and measurement updating, which is described below.

This starts with the following equations

$$z_{k|k-1}=z_{k-1|k-1}+{{\displaystyle \int }}_{t_{k-1}}^{t_k}{f}^c\left(z_{k-1|k-1},{\widehat{f}}_{k-1}\right)dt,$$

(8)

$${\widehat{f}}_k=J_k\left(y_k-h_k\left(z_{k|k-1},{\widehat{f}}_{k-1}\right)+{\nabla }_f{h}_k·{\widehat{f}}_{k-1}\right),$$

(9)

$$z_{k|k}=z_{k|k-1}+K_k\left[y_k-h\left(z_{k|k-1},{\widehat{f}}_{k-1}\right)-{\nabla }_f{h}_k·\left({\widehat{f}}_k-{\widehat{f}}_{k-1}\right)\right].$$

(10)

The matrices $J_k$ and $K_k$ need to be determined in order to estimate the state $z_{k|k}$ and the excitation ${\widehat{f}}_k$ as the minimum variance unbiased estimation of the system state $z_k$ and the real excitation $f_k$.

(1)

Step 1: State updating

The recursive relation of the extended state vector can be written as

$$\begin{array}{l}{z}_k=z_{k-1}+{{\displaystyle \int }}_{t_{k-1}}^{t_k}{f}^c\left(z\left(t\right),f\left(t\right)\right)dt+w=z_{k-1}+\Delta t·[f^c\left(z_{k-1|k-1},{\widehat{f}}_{k-1}\right)\\ +{\nabla }_z{f}_{k-1}^c·\left(z_{k-1}-z_{k-1|k-1}\right)+{\nabla }_f{f}_{k-1}^c·\left(f_{k-1}-{\widehat{f}}_{k-1}\right)]+w.\end{array}$$

(11)

Then, the prior estimation of the extended state vector can be written as

$$\begin{array}{l}{z}_{k|k-1}=z_{k-1|k-1}+{{\displaystyle \int }}_{t_{k-1}}^{t_k}{f}^c\left(z_{k-1|k-1},f_{k-1}\right)dt\\ \hspace{1em} =z_{k-1|k-1}+\Delta t·f^c\left(z_{k-1|k-1},f_{k-1}\right).\end{array}$$

(12)

Thus, the prior estimation error is

$$\begin{array}{l}{\tilde{z}}_{k|k-1}=z_k-z_{k|k-1}\\ \hspace{1em}\hspace{1em}=z_{k-1}-z_{k-1|k-1}+\Delta t·\left[{\nabla }_z{f}_{k-1}^c·\left(z_{k-1}-z_{k-1|k-1}\right)+{\nabla }_f{f}_{k-1}^c·\left(f_{k-1}-{\widehat{f}}_{k-1}\right)\right]+w\\ \hspace{1em}\hspace{1em}=\left(1+\Delta t·{\nabla }_z{f}_{k-1}^c\right)·\left(z_{k-1}-z_{k-1|k-1}\right)+\Delta t·{\nabla }_f{f}_{k-1}^c·\left(f_{k-1}-{\widehat{f}}_{k-1}\right)+w\\ \hspace{1em}\hspace{1em}=\left(1+\Delta t·{\nabla }_z{f}_{k-1}^c\right)·{\tilde{z}}_{k-1|k-1}+\Delta t·{\nabla }_f{f}_{k-1}^c·{\tilde{f}}_{k-1}+w.\end{array}$$

(13)

where ${\tilde{z}}_{k-1|k-1}$ is the posterior estimation error at time $(k-1)\Delta t$, and ${\tilde{f}}_{k-1}$ is the load estimation error at time $(k-1)\Delta t$.

We can obtain the variance of prior estimation error using

$$\begin{array}{l}{P}_{k|k-1}^z=E({\tilde{z}}_{k|k-1}·\left({\tilde{z}}_{k|k-1}{)}^T\right)\\ \hspace{1em}\hspace{1em}=\left[\begin{array}{cc}1+\Delta t·{\nabla }_z{f}_{k-1}^c& \Delta t·{\nabla }_f{f}_{k-1}^c\end{array}\right]\left[\begin{array}{cc}{P}_{k-1|k-1}^z& P_{k-1}^{zf}\\ P_{k-1}^{fz}& P_{k-1}^f\end{array}\right]\left[\begin{array}{c}{(1+\Delta t·{\nabla }_z{f}_{k-1}^c)}^T\\ {(\Delta t·{\nabla }_f{f}_{k-1}^c)}^T\end{array}\right]+Q.\end{array}$$

(14)

with $P_{k-1|k-1}^z$ being the variance matrix of a posterior error estimation at the time instant $(k-1)\Delta t$, while $P_{k-1}^f$ is the variance matrix of the load error at the same time instant. The values of $P_{k-1}^{zf}$ and $P_{k-1}^{fz}$ are the cross covariance matrix of the extended state vector and load:

$$P_{k-1|k-1}^z=E({\tilde{z}}_{k-1|k-1}·\left({\tilde{z}}_{k-1|k-1}{)}^T\right),$$

(15)

$$P_{k-1}^f=E({\tilde{f}}_{k-1}·\left({\tilde{f}}_{k-1}{)}^T\right),$$

(16)

$$P_{k-1}^{zf}={(P_{k-1}^{fz})}^T=E({\tilde{z}}_{k-1|k-1}·\left({\tilde{f}}_{k-1}{)}^T\right).$$

(17)

(2)

Step 2: Excitation identification

Firstly, we need to establish the equation about the load to be solved:

$${\tilde{y}}_k=y_k-h_k\left(z_{k|k-1},{\widehat{f}}_{k-1}\right)+{\nabla }_f{h}_k·{\widehat{f}}_{k-1}.$$

(18)

Then, by substituting Equation (7) into the above Equation (18), we get

$$\begin{array}{l}{\tilde{y}}_k={\nabla }_f{h}_k·f_k+{\nabla }_z{h}_k·\left(z_k-z_{k|k-1}\right)+v\\ \hspace{1em}={\nabla }_f{h}_k·f_k+{\nabla }_z{h}_k·{\tilde{z}}_{k|k-1}+v={\nabla }_f{h}_k·f_k+\tilde{e}\end{array}$$

(19)

where $\tilde{e}={\nabla }_z{h}_k·{\tilde{z}}_{k|k-1}+v$, and the mean value of $\tilde{e}$ is 0. As the components of $\tilde{e}$ have different variances, it is necessary to select the appropriate weighting matrix $W$ to optimize the estimation result. Then,

$${\widehat{f}}_k={[{({\nabla }_f{h}_k)}^T·W·{\nabla }_f{h}_k]}^{-1}{({\nabla }_f{h}_k)}^{T}W·{\tilde{y}}_k,$$

(20)

$${\tilde{R}}_k=E\left({\tilde{e}}_k{\tilde{e}}_k^T\right)={({\nabla }_z{h}_k)}^T·P_{k|k-1}^z·{\nabla }_z{h}_k+R_k.$$

(21)

where $W$ is the inverse matrix of the variance matrix ${\tilde{R}}_k$, and ${\tilde{R}}_k$ is the variance matrix of $\tilde{e}$. The matrix $J_k$ is

$$J_k={[{({\nabla }_f{h}_k)}^T{\tilde{R}}_k{}^{-1}{\nabla }_f{h}_k]}^{-1}{({\nabla }_f{h}_k)}^T{\tilde{R}}_k{}^{-1}.$$

(22)

In this case, the variance matrix of the error of the excitation estimation ${\widehat{f}}_k$ is as follows:

$$P_k^f={[{({\nabla }_f{h}_k)}^T·{\tilde{R}}_k{}^{-1}·{\nabla }_f{h}_k]}^{-1}.$$

(23)

(3)

Step 3: Measurement updating

After determining the matrix $J_k$, we also need to determine the Kalman gain matrix $K_k$. Substituting Equation (23) into Equation (11), we obtain the posterior estimation $z_{k|k}$:

$$\begin{array}{l}{z}_{k|k}=z_{k|k-1}+K_k[y_k-h\left(z_{k|k-1},{\widehat{f}}_{k-1}\right)\\ \hspace{1em} -{\nabla }_f{h}_k·\left(J_k\left(y_k-h_k\left(z_{k|k-1},{\widehat{f}}_{k-1}\right)+{\nabla }_f{h}_k·{\widehat{f}}_{k-1}\right)-{\widehat{f}}_{k-1}\right)]\\ \hspace{1em} =z_{k|k-1}+K_k\left(I-{\nabla }_f{h}_k·J_k\right)\left(y_k-h\left(z_{k|k-1},{\widehat{f}}_{k-1}\right)+{\nabla }_f{h}_k·{\widehat{f}}_{k-1}\right).\end{array}$$

(24)

Thus, the error of posterior estimation ${\tilde{z}}_{k|k}$ is

$$\begin{array}{l}{\tilde{z}}_{k|k}=z_k-z_{k|k}=z_k-z_{k|k-1}-K_k\left(I-{\nabla }_f{h}_k·J_k\right)\left(y_k-h\left(z_{k|k-1},{\widehat{f}}_{k-1}\right)+{\nabla }_f{h}_k·{\widehat{f}}_{k-1}\right)\\ =z_k-z_{k|k-1}-K_k\left(I-{\nabla }_f{h}_k·J_k\right)\left({\nabla }_z{h}_k·\left(z_k-z_{k|k-1}\right)+{\nabla }_f{h}_k·f_k+v\right)\\ =\left(I-K_k\left(I-{\nabla }_f{h}_k·J_k\right){\nabla }_z{h}_k\right)·{\tilde{z}}_{k|k-1}+K_k\left(I-{\nabla }_f{h}_k·J_k\right){\nabla }_f{h}_k·f_k+K_k\left(I-{\nabla }_f{h}_k·J_k\right)·v.\end{array}$$

(25)

Considering $z_{k|k-1}$ is an unbiased estimation, in order to ensure that $z_{k|k}$ is also an unbiased estimation, the following equation must be satisfied:

$$K_k\left(I-{\nabla }_f{h}_k·J_k\right){\nabla }_f{h}_k=0.$$

(26)

Therefore, the variance of posterior estimation can be written as

$$\begin{array}{l}{P}_{k|k}^z=E({\tilde{z}}_{k|k}·\left({\tilde{z}}_{k|k}{)}^T\right)\\ =\left[I-K_k\left(I-{\nabla }_f{h}_k·J_k\right){\nabla }_z{h}_k\right]·P_{k|k-1}^z·\left[I-K_k\left(I-{\nabla }_f{h}_k·J_k\right){\nabla }_z{h}_k{]}^T+K_k\left(I-{\nabla }_f{h}_k·J_k\right)·R·\right[K_k\left(I-{\nabla }_f{h}_k·J_k\right){]}^T\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =\left(I-L_k{\nabla }_z{h}_k\right)·P_{k|k-1}^z·{(I-L_k{\nabla }_z{h}_k)}^T+L_k·R·L_k{}^T.\end{array}$$

(27)

Among them,

$$L_k=K_k\left(I-{\nabla }_f{h}_k·J_k\right).$$

(28)

In order to make the posterior estimation optimal, the variance matrix $P_{k|k}^z$ of the posterior error estimation should be minimized:

$$\frac{\partial P_{k|k}^z}{\partial K_k}=0.$$

(29)

According to the above Equation (29), the Kalman filter gain $K_k$ can be obtained as follows:

$$K_k=P_{k|k-1}^z{({\nabla }_z{h}_k)}^T{\tilde{R}}_k{}^{-1}.$$

(30)

By substituting the above Equation (30) into Equation (27), the variance matrix of posterior error estimation is rewritten as

$$P_{k|k}^z=P_{k|k-1}^z-K_k·({\tilde{R}}_k-{\nabla }_f{h}_k·P_k^f·\left({\nabla }_f{h}_k{)}^T\right)·K_k.$$

(31)

Then, we can get

$$P_k^{zf}=-P_{k|k-1}^z{({\nabla }_z{h}_k)}^T{J}_k{}^T=-K_k·{\nabla }_f{h}_k·P_k^f.$$

(32)

In order to reduce the numbers of the given initial conditions, the order of the load identification algorithm based on the extended Kalman filter can be adjusted including all the above three steps. It should be stressed that the order adjustment does not affect the efficiency of the algorithm. It only involves a different initial value which has little influence on the whole process of load identification.

In this paper, we distinguish between two cases:

• the unknown parameters are the stiffness or the damping coefficients
• the unknown parameter are the mass coefficients.

First, we discuss the case of unknown stiffness or damping. We can get the partial derivative matrix of the state update equation and the measurement update equation to the extended state vector based on Equations (4) and (5) such as

$$\begin{array}{c}{\nabla }_z{f}_{k-1}^c=\left[\frac{\partial f^c}{\partial p} \frac{\partial f^c}{\partial \dot{p}} \frac{\partial f^c}{\partial \alpha }\right]\left|z=z_{k-1\left|k-1\right.}\right.\\ =M^{-1}\left[\begin{array}{ccccc}\mathbf{0}& \mathit{I}& \mathbf{0}& \cdots & \mathbf{0}\\ -K& -C& -\frac{\partial K}{\partial {\alpha }_1}{p}_{k-1\left|k\right.-1}-\frac{\partial C}{\partial {\alpha }_1}{\dot{p}}_{k-1\left|k-1\right.}& \cdots & -\frac{\partial K}{\partial {\alpha }_m}{p}_{k-1\left|k-1\right.}-\frac{\partial C}{\partial {\alpha }_m}{\dot{p}}_{k-1\left|k-1\right.}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \cdots & \mathbf{0}\end{array}\right] \end{array}$$

(33)

$$\begin{array}{c}{\nabla }_z{h}_k=\left[\begin{array}{ccc}\frac{\partial h}{\partial p}& \frac{\partial h}{\partial \dot{p}}& \frac{\partial h}{\partial \alpha }\end{array}\right]\left|z=z_{k\left|k-1\right.}\right.\\ =M^{-1}{H}_0\left[\begin{array}{ccc}-K& -C& -\frac{\partial K}{\partial {\alpha }_1}{p}_{k\left|k-1\right.}-\frac{\partial C}{\partial {\alpha }_1}{\dot{p}}_{k\left|k-1\right.}\end{array}\left.\begin{array}{cc}\cdots & -\frac{\partial K}{\partial {\alpha }_s}{p}_{k\left|k-1\right.}-\frac{\partial C}{\partial {\alpha }_s}{\dot{p}}_{k\left|k-1\right.}\end{array}\right]\right.\end{array}$$

(34)

The partial derivative matrices of the state update equation and the measurement update equation with respect to the load vector are

$${\nabla }_f{f}_{k-1}^c=\left[\begin{array}{c}\mathbf{0}\\ \begin{array}{c}{M}^{-1}{B}_u\\ \mathbf{0}\end{array}\end{array}\right],$$

(35)

$${\nabla }_f{h}_k=H_0{M}^{-1}{B}_u.$$

(36)

In this case, ${\nabla }_f{f}_{k-1}^c \mathrm{and} {\nabla }_f{h}_k$ are constant matrices. The measurement update step and the load identification step can be simplified as

$$z_{k|k}=z_{k|k-1}+K_k\left[y_k-\bar{h}\left(z_{k|k-1}\right)-{\nabla }_f{h}_k·{\widehat{f}}_k\right],$$

(37)

$${\widehat{f}}_k=J_k\left(y_k-{\bar{h}}_k\left(z_{k|k-1}\right)\right).$$

(38)

Thus, the algorithm flow of dynamic load identification can be obtained.

Second is the case where the mass coefficients are unknown. We define the symmetric matrix $M^{-1}$ as

$$M^{-1}=\left[\begin{array}{cccc}\frac{1}{m_1}& & & \\ & \frac{1}{m_2}& & \\ & & \ddots & \\ & & & \frac{1}{m_n}\end{array}\right].$$

(39)

All n mass coefficients of the structure are then inserted in the extended state vector. We can thus get the partial derivatives of the state updating equation and the measurement updating equation as

$${\nabla }_z{f}_{k-1}^c=\left[\begin{array}{ccc}\frac{\partial f^c}{\partial p}& \frac{\partial f^c}{\partial \dot{p}}& \frac{\partial f^c}{\partial \alpha }\end{array}\right]\left|z=z_{k\left|k\right.}\right.\hspace{0.33em}=\left[\begin{array}{ccc}\mathbf{0}& I& \mathbf{0}\\ -M^{-1}K& -M^{-1}C& E_k\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\end{array}\right],$$

(40)

$${\nabla }_z{h}_k=\left[\begin{array}{ccc}\frac{\partial h}{\partial p}& \frac{\partial h}{\partial \dot{p}}& \frac{\partial h}{\partial \alpha }\end{array}\right]\left|z=z_{k\left|k-1\right.}\right.\hspace{0.33em}=\left[\begin{array}{ccc}-H_0{M}^{-1}K& -H_0{M}^{-1}C& H_0{\bar{E}}_k\end{array}\right].$$

(41)

The symmetric matrices $E_k$ and ${\bar{E}}_k$ can be written as

$$\begin{array}{l}{E}_k=\left[\begin{array}{cccc}-\frac{1}{m_1^2}{\xi }_k\left(1\right)& & & \\ & -\frac{1}{m_2^2}{\xi }_k\left(2\right)& & \\ & & \ddots & \\ & & & -\frac{1}{m_n^2}{\xi }_k\left(n\right)\end{array}\right],\\ {\bar{E}}_k=\left[\begin{array}{cccc}-\frac{1}{m_1^2}{\bar{\xi }}_k\left(1\right)& & & \\ & -\frac{1}{m_2^2}{\bar{\xi }}_k\left(2\right)& & \\ & & \ddots & \\ & & & -\frac{1}{m_n^2}{\bar{\xi }}_k\left(n\right)\end{array}\right].\end{array}$$

(42)

Let ${\xi }_k(i)$, ${\bar{\xi }}_k(i)$ be the i-th element of the vector ${\xi }_k$, ${\bar{\xi }}_k$ at time $(k-1)\Delta t$

$${\xi }_k=-K{p}_{k|k}-C{p}_{k|k}+B_u{f}_k,$$

(43)

$${\bar{\xi }}_k=-K{p}_{k|k-1}-C{p}_{k|k-1}+B_u{f}_{k-1}.$$

(44)

The partial derivatives of the state updating equation and the measurement updating equation with respect to the load vector are the same as the Equations (35) and (36). However, they are not constant matrices in this case because the mass coefficient will change in the filtering process. The mass coefficient must be replaced back to the mass matrix in each filtering. The algorithm procedure is then the same as that of the unknown parameter stiffness or damping coefficient, but the difference is in calculating the Jacobian matrix.

For simplicity, we summarize the algorithm steps of the load identification in Algorithm 1.

Algorithm 1. Flow of load identification based on extended Kalman filter

1. Given initial conditions               $z_{0|-1},P_{0|-1}^{z}and f_0$   2. Excitation identification step            ${\tilde{R}}_k={({\nabla }_z{h}_k)}^T·P_{k|k-1}^z·{\nabla }_z{h}_k+R$           $J_k={[{({\nabla }_f{h}_k)}^T{\tilde{R}}_k{}^{-1}{\nabla }_f{h}_k]}^{-1}{({\nabla }_f{h}_k)}^T{\tilde{R}}_k{}^{-1}$          ${\widehat{f}}_k=J_k\left(y_k-h_k\left(z_{k|k-1},{\widehat{f}}_{k-1}\right)+{\nabla }_f{h}_k·{\widehat{f}}_{k-1}\right)$             $P_k^f={[{({\nabla }_f{h}_k)}^T·{\tilde{R}}_k{}^{-1}·{\nabla }_f{h}_k]}^{-1}$   3. Measurement update step             $K_k=P_{k|k-1}^z{({\nabla }_z{h}_k)}^T{\tilde{R}}_k{}^{-1}$      $z_{k|k}=z_{k|k-1}+K_k\left[y_k-h\left(z_{k|k-1},{\widehat{f}}_{k-1}\right)-{\nabla }_f{h}_k·\left({\widehat{f}}_k-{\widehat{f}}_{k-1}\right)\right]$        $P_{k|k}^z=P_{k|k-1}^z-K_k·({\tilde{R}}_k-{\nabla }_f{h}_k·P_k^f·\left({\nabla }_f{h}_k{)}^T\right)·K_k$            $P_{k|k}^{zf}={(P_{k|k}^{fz})}^T=-K_k·{\nabla }_f{h}_k·P_k^f$   4. State update step              $z_{k+1|k}=z_{k|k}+{{\displaystyle \int }}_{t_k}^{t_{k+1}}{f}^c\left(z_{k|k},{\widehat{f}}_k\right)dt$      $P_{k+1|k}^z=\left[\begin{array}{cc}1+\Delta t·{\nabla }_z{f}_k^c& \Delta t·{\nabla }_f{f}_k^c\end{array}\right]\left[\begin{array}{cc}{P}_{k|k}^z& P_k^{zf}\\ P_k^{fz}& P_k^f\end{array}\right]\left[\begin{array}{c}{(1+\Delta t·{\nabla }_z{f}_k^c)}^T\\ {(\Delta t·{\nabla }_f{f}_k^c)}^T\end{array}\right]+Q$

3. Model Reduction Strategy

The proposed algorithm works online to recover the applied load and the system characteristics as they change in real-time, However, when considering a large system with thousands of degrees of freedom, it becomes impossible to perform the computations in real-time. In this case, the model reduction technique can be used to reduce the involved computations. Here, we adapt the reduction strategy. The main idea in this strategy is to divide the degrees of freedom of the finite element model into master and slave degrees of freedom. The master degrees of freedom are generally those that have forces applied at them, or we need to read the results at them. The other degrees of freedom can be defined as slave degrees of freedom. The slave degrees of freedom are expressed in terms of the master degree of freedom which significantly reduces the required computations.

Next, we explain the reduction strategy using a beam finite element with four degrees of freedom associated to two nodes so that each node has a rotational and a translational degree of freedom. The rotational degrees of freedom are difficult to measure, while the translational ones are easy. Using the reduction method, the dynamic equations of the beam elements can be transformed into condensed dynamic equations with only translational degrees of freedom.

The mass associated to the translational degrees of freedom is considered and is assumed to be concentrated, while the rotational degrees of freedom mass are ignored. In this case, the stiffness matrix and the mass matrix of the beam element can be rewritten as

$$\begin{array}{llll}& u_1{u}_2{\theta }_1{\theta }_2& & u_1{u}_2{\theta }_1{\theta }_2\\ K_i=\frac{E_i{I}_i}{l^3}·& \left[\begin{array}{cccc}\hspace{0.33em}12& -12& 6l& \hspace{0.33em}\hspace{0.33em}6l\\ -12& 12& -6l& -6l\\ 6l& -6l& \hspace{0.33em}\hspace{0.33em}4l^2& \hspace{0.33em}\hspace{0.33em}2l^2\\ 6l& -6l& \hspace{0.33em}\hspace{0.33em}2l^2& \hspace{0.33em}\hspace{0.33em}4l^2\end{array}\right]& \hspace{0.33em}, \hspace{0.33em}{M}_i=\frac{{\rho }_i{A}_{i}l}{2}·& \left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right].\end{array}$$

(45)

The global stiffness and mass matrices are composed of their counterpart elementary matrices. The first half of the global stiffness and mass matrices correspond to the translational degrees of freedom, while the second half corresponds to the rotational degrees of freedom. Then, the global stiffness matrix and mass matrix can be written as

$$K=\left[\begin{array}{cc}{K}_{uu}& K_{u\theta }\\ K_{\theta u}& K_{\theta \theta }\end{array}\right], M=\left[\begin{array}{cc}{M}_{uu}& 0\\ 0& 0\end{array}\right].$$

(46)

Substituting the global stiffness and mass matrices into the differential equation of motion, we can get

$$\left[\begin{array}{cc}{M}_{uu}& 0\\ 0& 0\end{array}\right]\left(\begin{array}{c}{\ddot{X}}_u\\ {\ddot{X}}_{\theta }\end{array}\right)+\left\{\alpha \left[\begin{array}{cc}{M}_{uu}& 0\\ 0& 0\end{array}\right]+\beta \left[\begin{array}{cc}{K}_{uu}& K_{u\theta }\\ K_{\theta u}& K_{\theta \theta }\end{array}\right]\right\}\left(\begin{array}{c}{\dot{X}}_u\\ {\dot{X}}_{\theta }\end{array}\right)+\left[\begin{array}{cc}{K}_{uu}& K_{u\theta }\\ K_{\theta u}& K_{\theta \theta }\end{array}\right]\left(\begin{array}{c}{X}_u\\ X_{\theta }\end{array}\right)=\left(\begin{array}{c}F\left(t\right)\\ 0\end{array}\right).$$

(47)

Among them, the parameter $X_u$ and $X_{\theta }$ mean translational displacement and rotational displacement, respectively. We also assume that the external load only acts on the beam structure in the translational direction. Thereafter, Equation (47) can be divided into two equations:

$$M_{uu}{\ddot{X}}_u+\left(\alpha M_{uu}+\beta K_{uu}\right){\dot{X}}_u+\beta K_{u\theta }{\dot{X}}_{\theta }+K_{uu}{X}_u+K_{u\theta }{X}_{\theta }=F\left(t\right),$$

(48)

$$\beta K_{\theta u}{\dot{X}}_u+\beta K_{\theta \theta }{\dot{X}}_{\theta }+K_{\theta u}{X}_u+K_{\theta \theta }{X}_{\theta }=0.$$

(49)

We can get

$$\beta {\dot{X}}_{\theta }+X_{\theta }=-K_{\theta \theta }^{-1}{K}_{\theta u}\left(\beta {\dot{X}}_u+X_u\right).$$

(50)

Substituting Equation (50) into Equation (48), we can get the motion equation of reduction

$$M_{gu}·{\ddot{X}}_u{+\mathrm{C}}_{gu}·{\dot{X}}_u+K_{gu}·X_u=F\left(t\right).$$

(51)

where

$$M_{gu}=M_{uu},$$

(52)

$$C_{gu}=\alpha M_{uu}+\beta \left(K_{uu}-K_{u\theta }{K}_{\theta \theta }^{-1}{K}_{\theta u}\right),$$

(53)

$$K_{gu}=K_{uu}-K_{u\theta }{K}_{\theta \theta }^{-1}{K}_{\theta u}.$$

(54)

Substitute the condensed dynamic equations into the proposed dynamic load identification algorithm. It can be seen that the process of load identification is the same as that without the reduction. Similar to before and after the reduction, the state update equation and the measurement update equation in the algorithm can be rewritten as

$$\frac{d}{dt}\left[\begin{array}{c}{X}_u\\ \begin{array}{c}{\dot{X}}_u\\ \alpha \end{array}\end{array}\right]=\left[\begin{array}{c}{\dot{X}}_u\\ M_{gu}{}^{-1}\left(-K_{gu}{X}_u-C_{gu}{\dot{X}}_u+B_{u}f\left(t\right)\right)\\ \mathbf{0}\end{array}\right]\hspace{0.33em}=f^c\left(z\left(t\right),f\left(t\right)\right),$$

(55)

$$y\left(t\right)=H_0{M}_{gu}{}^{-1}\left(-C_{gu}{\dot{X}}_u-K_{gu}{X}_u+B_{u}f\left(t\right)\right)=h\left(z\left(t\right),f\left(t\right)\right).$$

(56)

4. Numerical Validation

In this section, three numerical examples are used to validate the proposed algorithm. The first two represent systems of multi-degrees of freedom while the third a beam. Two of the examples are dedicated to systems with unknown mass while the third to unknown stiffness. The measurement taken in the examples is assumed to be contaminated with low and high levels of noise. In the last example, the proposed algorithm is combined with the reduction algorithm to identify the system and its load.

4.1. Three-Degrees-of-Freedom with Varying Mass

In the first example, a three-degrees-of-freedom system is selected. A schematic diagram of the structure is shown in Figure 1. The damping of the system is assumed to be Rayleigh damping, with the damping coefficients $\alpha$ = 0.05, $\beta$ = 0.02.

Here, m₁ = m₂ = m₃ = 1 kg, k₁ = k₂ = k₃ = k₄ = 200 N/m. We assume the mass coefficient m₂ increases gradually from 1 kg to 3 kg starting 1.5 s and finishing at 3.5 s. The responses of the three masses are selected as the measurement vector, while the three masses are included in the extended state vector as unknown parameters. Assuming that the estimation of the initial response is correct and the estimation of the initial parameters is biased, the initial extended state vector can be set as $z_{0|-1}={[0,0,0,0,0,0,1,3,4]}^T$, which is a 9 × 1 column vector. We set the covariance matrix of the model noise to $Q=diag\left(\left(1\times {10}^{-6}\mathrm{ones}\left(1,6\right),1\times {10}^{-4}\mathrm{ones}\left(1,3\right)\right)\right),$ and the covariance matrix of the measurement noise to $R=1\times {10}^{-6}\mathrm{eye}\left(3\right)$.

A load of $f=sin(5\pi t)+2sin(2\pi t)$ is applied on the mass block m₁ and 5%, 10% noise are applied to the measurement vector, respectively. The identification results are shown in the Figure 2, Figure 3, Figure 4 and Figure 5. The left plots are for the case with 5% noise, while the right side is for the case with 10% noise.

The results show that the proposed algorithm can accurately and simultaneously identify the varying masses and the unknown loads on the structure when the noise level is 5% and 10%. The values of the unknown parameters identified by the algorithm are recorded before the mass changes (t = 1 s) and after it changes (t = 5 s), and are compared with the real values in

Table 1. Parameter identification error.

Mass Coefficient		t = 1 s			t = 5 s
		Identification (kg)	Exact (kg)	Error (%)	Identification (kg)	Exact (kg)	Error (%)
5% noise	m1	1.013	1	1.3	0.9867	1	1.3
	m2	1.087	1	8.7	3.022	3	0.7
	m3	0.9856	1	1.4	1.006	1	0.6
10% noise	m1	1.002	1	0.2	1.142	1	14.2
	m2	1.025	1	2.5	3.049	3	4.9
	m3	1.051	1	5.1	1.139	1	13.9

. The identification errors are shown in the same table. It can be seen that the algorithm has successfully identified the unknown mass values.

In order to illustrate the identification accuracy of the unknown loads, the relative error (RE) and the correlation coefficient (r) methods are used to verify the results. The relative error refers to the ratio of the absolute error between the identification result and the exact value divided by the exact value. The smaller the relative error is, the higher the accuracy of the identification result. The correlation coefficient (r) is the quantity to study the level of the linear correlation between the variables and is expressed as a percentage. The closer the correlation coefficient is to 1, the better the identification. The errors are shown in

Table 2. Load identification error.

Noise Level	Relative Error (RE) (%)	Correlation Coefficient (r) (%)
5% noise	9.55	99.48
10% noise	17.75	98.44

. It can be seen that the algorithm can accurately identify the unknown load in general even when the measured data is contaminated with 5% or 10% noise.

4.2. Five-Degree-of-Freedom with Varying Stiffness

In the second numerical test, we study a five-degrees-of-freedom system. We show in Figure 6 a schematic diagram of the structure. In this test, we consider a gradually varying stiffness. The damping of the system is assumed to be Rayleigh damping of the coefficients $\alpha$ = 0.05, $\beta$ = 0.02.

The considered masses and stiffnesses are m₁ = m₂ = m₃ = m₄ = m₅ = m₆ = 1 kg and k₁ = k₂ = k₃ = k₄ = k₅ = k₆ = 200 N/m. The stiffness coefficient k₄ of the structure decreases gradually from 200 N/m to 120 N/m starting at the time instant t = 1.5 s and finishes at t = 3.5 s. The response of all the mass blocks is recorded in the measurement vector. Assuming that the stiffness coefficients k₁ and k₂ are accurately estimated, the remaining stiffness coefficients are substituted into the extended state vector as unknown parameters. We start by considering the initial response estimation to be correct and that the initial parameter estimation has a certain deviation. The initial extended state vector can be set to $z_{0|-1}={(zeros\left(1,10\right),\hspace{0.33em}120,220,160,180)}^T$, which is a 14 × 1 column vector.

The mass m₁ is excited with a load given by $f_1=sin(5\pi t)+2sin(2\pi t),$ while the mass m₂ by a load $f_2=sin(6\pi t)$. We consider the measurement to be polluted with noise so that we first add 1% and then 5% white noise to the measurement vector. Again, we use the proposed algorithm to track the unknown parameters as well as identify the unknown load. Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11 show the identified stiffness while Figure 12 and Figure 13 show the identified loads. The left plots show the identification results when the measurement is contaminated with 1% noise, while the right plots with 5% noise. The values of the unknown parameters identified by the algorithm are recorded before the parameters change at t = 1 s and again after the parameters change at t = 5 s, and the results are compared to the real values in

Table 3. Parameter identification error.

Stiffness Coefficient		t = 1 s			t = 5 s
		Identification (N/m)	Exact (N/m)	Error (%)	Identification (N/m)	Exact (N/m)	Error (%)
1% noise	k3	198.8	200	0.6	199.3	200	0.35
	k4	199.2	200	0.4	118.4	120	1.33
	k5	199.6	200	0.2	198.9	200	0.55
	k6	199.8	200	0.1	199.2	200	0.4
5% noise	k3	187.2	200	6.4	194.4	200	2.8
	k4	191.8	200	4.1	118.8	120	1
	k5	198.7	200	0.65	197.5	200	1.25
	k6	201	200	0.5	195.6	200	2.2

. The corresponding errors are also shown in the table. The errors in the identified load under different noise conditions are shown in the

Table 4. Load identification error.

Noise Level		Relative Error (RE) (%)	Correlation Coefficient (r) (%)
1% noise	f1	2.64	99.97
	f2	4.98	99.88
5% noise	f1	10.24	99.32
	f2	19.85	98.05

.

The results show that the proposed algorithm can successfully identify the unknown load with varying stiffness even when the measurement is polluted with noise. The algorithm is also capable of tracking the changes in the stiffness as it recovers the load. The results suggest that the identification results of the proposed algorithm are accurate.

4.3. Cantilever with Varying Mass

In the third example, we study a cantilever with a varying mass. The considered beam length is L = 1 m with a rectangular cross section of the width b = 0.1 m and the height h = 0.01 m. The material density is $\rho =2770 \mathrm{kg}/{\mathrm{m}}^3$. The Young’s modulus is E = 71 GPa. The damping characteristic of the cantilever is assumed to be Rayleigh damping with the coefficients being $\alpha =0.8$, $\beta =0.004$. We mesh the cantilever using six beam finite elements. Figure 13 shows the cantilever and the considered mesh. In this example, we consider the mass ${\rho }_i{A}_i$ of the beam element to be gradually changing. We assume the mass characteristics ${\rho }_4{A}_4$ of the element number 4 to decrease linearly between 1.5 s and 3.5 s to 60% of its original value. The masses of the other elements do not change. A sinusoidal force $F=200sin(5\pi t)$ is applied at the end of the cantilever. Again, two cases are considered where 1% and then 5% white Gaussian noise is added to the measured response.

The reduction method proposed in Section 3 is used to reduce the dynamic equation of the beam structure to a dynamic equation with only translational degrees of freedom. Then, the proposed algorithm is used to identify the load and the parameters of the condensed beam structure. We assume that the mass coefficient ${\rho }_6{A}_6$ is accurately estimated, while other elements’ mass characteristics are assumed unknown. Thus, the extended state vector of the proposed algorithm can be written as $z={(p,\dot{p},{\rho }_1{A}_1,{\rho }_2{A}_2,{\rho }_3{A}_3,{\rho }_4{A}_4,{\rho }_5{A}_5)}^T$. Let the initial extended state vector be $z_{0|-1}={(zeros\left(1,12\right),\hspace{0.33em}\rho A,0.6\rho A.\mathrm{ones}\left(1,4\right))}^T$, which is a 17 × 1 column vector, and all the translational degrees of freedom on the beam are observed. Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18 show the identification results of the unknown mass parameters, while Figure 19 shows the identified load. Again, on the left side, we show the results when 1% noise is applied to the measurement, while to the right are the results correspond to 5% noise pollution.

The results’ accuracy clearly shows that the reduction algorithm works well with the proposed load identification algorithm. The load as well as the unknown mass of each beam element are recovered with good accuracy. The errors of the identification results are shown in

Table 5. Parameter identification error.

Mass Characteristic		t = 1 s			t = 5 s
		Identification (kg/m)	Exact (kg/m)	Error (%)	Identification (kg/m)	Exact (kg/m)	Error (%)
1% noise	${\rho }_1{A}_1$	2.76	2.77	3.61	2.903	2.77	4.8
	${\rho }_2{A}_2$	2.818	2.77	1.73	2.506	2.77	9.53
	${\rho }_3{A}_3$	2.765	2.77	0.18	2.944	2.77	6.28
	${\rho }_4{A}_4$	2.773	2.77	0.11	1.554	1.662	6.49
	${\rho }_5{A}_5$	2.785	2.77	0.54	2.908	2.77	4.98
5% noise	${\rho }_1{A}_1$	2.773	2.77	0.11	2.829	2.77	2.13
	${\rho }_2{A}_2$	2.285	2.77	17.51	2.349	2.77	15.20
	${\rho }_3{A}_3$	3.172	2.77	14.51	3.155	2.77	13.89
	${\rho }_4{A}_4$	2.518	2.77	9.09	1.863	1.662	7.26
	${\rho }_5{A}_5$	2.924	2.77	5.56	2.949	2.77	6.46

and

Table 6. Load identification error.

Noise Level	Relative Error (RE) (%)	Correlation Coefficient (r) (%)
1% noise	2.74	99.96
5% noise	13.66	99.08

for the mass and the load, respectively. The results with both levels of noise are accurate with engineering accuracy, i.e., the error is the same level as that of the pollution noise. Obviously, increasing the noise level results in larger errors, but the algorithm remains to provide meaningful results.

5. Experimental Validation

As a final validation, an online dynamic load identification experiment of a simply supported beam with varying mass characteristics is designed to verify the proposed algorithm experimentally. The structural parameters of the simply supported beam are shown in

Table 7. Parameter table of simply supported beam.

Parameters of Simply Supported Beam	Value
Length $a$	0.7 m
Width $b$	0.04 m
Thickness $h$	0.008 m
Density $\rho$	7800 kg/m3
Elastic Modulus $E$	209 GPa
Poisson’s ratio $\mu$	0.30

.

Before we run the experiment, the natural frequencies of the simply supported beam are evaluated both experimentally and numerically using the finite element method. The natural frequencies and natural mode shapes of the beam are shown in

Table 8. Natural frequency comparison.

Order	Experimental Results (Hz)	Simulation Results (Hz)	Error (%)
1	37.864	37.8051	0.156
2	151.74	151.264	0.314
3	337.58	340.8637	0.963
4	603.5	609.2473	0.943

where the numerical and experimental results are compared. It can be seen that the natural frequencies of the experimental results are very close to that of the simulation model. This is reflected in the first four mode shapes where the relative errors are less than 1%. Thus, we conclude that the finite element model is reliable.

To vary the beam mass, we fix an aluminum box at 0.48 m from the left of the simply supported beam. The aluminum box is filled with steel sand. To change the beam mass, the box has holes at the bottom so that it leaks the steel sand as the experiment progresses. The aluminum box with the iron sand and the holes is shown in Figure 20, while a schematic diagram of the experiment is shown in Figure 21. The experimental setting is shown in Figure 22. All the used instruments are listed in

Table 9. Force identification experiment equipment.

Equipment Classification	Name
Vibration Exciter	JZT-2 Permanent magnet exciter
Power amplifier	HEAS-50 Power amplifier
Dynamic signal acquisition board	NI PXIe-4499 Capture card
Signal acquisition instrument	NI PXI
Sensor	PCB 356A33 accelerometer
Software	NI Signal Express

. The sampling rate is set to 8196 Hz, and the excitation signal on the exciter is measured by a force sensor. The beam is meshed into eight beam elements. The aluminum box is fixed to element 6. The acceleration responses are measured by sensors placed at the nodes as shown in Figure 21. Two loading cases are considered. First, a sinusoidal load is applied, and then the experiment is repeated for a triangular wave load.

Using the acceleration measurements at the nodes, the proposed identification algorithm is used to identify the unknown load and the unknown changing mass. The identification results are evaluated by comparing them to the actual load measured by the force sensor. The mass identification results and evaluation results are shown in Figure 23, Figure 24, Figure 25 and Figure 26.

Table 10. Parameter identification errors.

Mass Characteristic	Before Mass Changes			After Mass Changed
	Identification (kg/m)	Exact (kg/m)	Error (%)	Identification (kg/m)	Exact (kg/m)	Error (%)
${\rho }_1{A}_1$	2.097	2.528	17.05	2.168	2.528	14.24
${\rho }_2{A}_2$	2.395	2.528	5.26	2.540	2.528	0.47
${\rho }_3{A}_3$	2.246	2.528	11.15	2.451	2.528	3.04
${\rho }_4{A}_4$	2.557	2.528	1.14	2.477	2.528	2.02
${\rho }_5{A}_5$	2.596	2.528	2.69	2.622	2.528	3.72
${\rho }_6{A}_6$	5.582	5.385	3.66	2.728	2.528	7.91
${\rho }_7{A}_7$	2.588	2.528	2.37	2.562	2.528	1.34
${\rho }_8{A}_8$	2.457	2.528	2.81	2.503	2.528	0.98

shows the error in the identified mass of individual elements. The errors are shown before the box has started leaking the sand and then again immediately after all the sand has leaked.

It can be seen that the proposed algorithm can identify the unknown mass with good accuracy regarding the parameters of the structure as well as the unknown loads acting on it. Other than the identification error of the mass characteristic ${\rho }_1{A}_1$, the identification of other mass coefficients is good. The errors between the identification results and the real values are due to the noise in the measurement responses and the inaccuracy caused by the approximation of the reduction method. Since most of the identification errors of these mass characteristic is within a close range of the element mass, it can be considered that the parameters identified by the proposed algorithm are accurate. Figure 27 and Figure 28 show the identification results of the sinusoidal loads before and during the beam mass change. While Figure 29 and Figure 30 show the identification results of the triangular loads before and during the beam mass change. In general, it can be seen that the algorithm has captured the loads with good accuracy. To quantify the identification accuracy further, we show in

Table 11. Load identification errors.

Type	Relative Error (RE) (%)	Correlation Coefficient (r) (%)
sine	9.31	99.57
triangular wave	12.97	99.16

the relative errors and the correlation coefficient for the identified results against the experimental results.

It can be seen in the figures that the accuracy of the load identification is indeed accurate where all the details of the applied load are captured. However, in the initial stage for the first 0.05 s, there are fluctuations in the results. This is caused by the initial errors in estimating the mass coefficients at this initial phase. With the identification progress, the mass parameter converges to accurate values, and the load identification accuracy significantly improves. The errors in Table 11 show the correlation coefficients of the sinusoidal load and the triangular load to be more than 99%, which suggests an excellent overall accuracy in the algorithm. This accuracy in identifying the dynamic load and in capturing the mass parameters has been achieved despite the varying nature of the problem.

6. Conclusions

In this paper, an online dynamic load identification method based on an extended Kalman filter is proposed. The proposed algorithm is also combined with the reduction method which is used to reduce the number of degrees of freedom in the system. The reduced system leads to improved efficiency by reducing the rank of the matrix to be solved, hence, the amount of necessary computation. The algorithm is implemented for cases when the unknown parameters are either the stiffness or the mass of a structure. It can be seen that the different unknown parameter types only affect the Jacobian matrix, but the algorithm flow in a similar way for both cases.

Several examples are used to test and validate the proposed approach. The studied examples include multi-degrees of freedom systems as well as a cantilever where the measured data are contaminated with different levels of noise. The algorithm was also validated in a lab experiment designed to validate the approach in a physical test on a simply supported beam with a varying mass. The reduction method is successfully implemented with the proposed algorithm to eliminate the rotational degrees of freedom in the beam elements. The numerical as well as the experimental results prove that the proposed algorithm can simultaneously identify the change in the structural parameters and capture them as well as identify the unknown dynamic load acting on the structure. All the presented results show the accuracy of the proposed algorithm and its efficiency in capturing the dynamic loads as well as the changing parameters in the system even when the measurement data are polluted with noise. Finally, for future works, we aim to develop a theoretical proof for the method converge. However, developing such a theoretical proof can be especially difficult given the varying nature of the identified parameters.