Real-Time Identification of Dynamic Loads Using Inverse Solution and Kalman Filter

Evaluating dynamic loads in real time is crucial for health monitoring, fault diagnosis and fatigue analysis in aerospace, automotive and earthquake engineering among other vibration related applications. Developing such algorithms can be vital for several safety and performance functionalities. Therefore, over the past few years the identification of dynamic loads has attracted a lot of attention; however, little literature on the online identification can be found. In this paper, we propose an online-identification method of structural dynamic loads so that the dynamic load is evaluated in real time and while the system response is still being measured. This is achieved by significantly improving the identification efficiency while retaining a high accuracy. The proposed method which is based on Kalman filter, is introduced in detail for a finite as well as an infinite number of degrees of freedom. Starting from an initial guess of the state vector we evaluate the error covariance, which then helps to identify the value of the excitation force using a weighted least square method and minimizing the covariance unbiased estimation. This is repeated at certain time intervals i.e., time steps where the state vector is updated in real time as acceleration measurements are updated. The feasibility of the method is validated using numerical simulations and an experimental verification where a detailed LabVIEW (National Instruments Ltd.) implementation is provided.


Introduction
Dynamic load identification is a relatively wide area of research that covers a range of applications [1][2][3][4][5]. It is crucial to evaluate the dynamic loads on a mechanical system in many applications such as design optimization, diagnostics, control and monitoring. It is common to record the system response in terms of displacement, velocity or acceleration which can then be processed offline to evaluate the external forces/loads. Although this is a relatively well-established area of research but it is still under continuous development where several challenges remains to be unsolved [2,6,7]. A major challenge in this regards is about moving the load identification process from offline to online computations so that the dynamic load is evaluated in real time and while the system response is still being measured.
In this paper, we focus on identifying the dynamic load in real time which we refer to as the online load identification. This challenge is mainly motivated by the automotive industry and is highly relevant for applications in autonomous vehicles [8][9][10]. For example in autonomous vehicles it is important to evaluate the external dynamic load on a vehicle to take decisions related to acceleration/deceleration or trajectory optimization [9]. Having real-time accurate estimates of the dynamic loads on such systems can help to prevent catastrophic failures. also considered in space with Kalman filter [42,43]. Although there are some advantages for these identification methods based on Kalman filter in some specific cases, compared to our proposed method, these different methods were only adopted for the discrete systems with few degrees of freedom, without considering the running of the loads in real application, and cannot reconstruct the external load in real time.
Different load identification methods were also developed based on the idea of Kalman filter. In order to identify the external excitation acting on a non-linear structural system the extended Kalman filter was combined with the recursive least square method [44] or the unscented Kalman filter is used [38]. The latter algorithm was verified by comparing the results to those found using numerical simulations of a non-linear centralized mass system as well as using experimental results. Moreover, the filter was also used for identifying excitation with a linear time discretization based on the minimum-variance unbiased estimate of the least square method [45]. The proposed algorithm is globally optimal in the minimum-variance unbiased sense. In a different application Kalman filter was used to separate the modal wind load from the measured acceleration response [46]. The modal wind load on a structure were identified by establishing the relevant filter gain matrix and using the noise covariance matrix. Moreover, for excitation identification on a non-linear structure, a method combining the least squares algorithm with adaptive decay factor was proposed [47] where a dual Kalman filter method was developed for identifying excitations. Here a numerical model is used to verify the anti-drift performance of the method.
In this study, we aim to estimate a dynamic load in the time domain by introducing an excitation identification step in Kalman filter. This step combines the time update and the measurement update steps to recover the external force in real time. The proposed method works with a finite as well as an infinite number of degrees of freedom using the modal transform method. The method is validated using experimental results as well as by comparison with a standard method. The source code of the developed algorithm is implemented in LabVIEW (National Instruments Ltd.) and made available in the paper Appendix A. The implementation is used to verify the algorithm experimentally. The main objective of our work is to provide an online dynamic load identification approach by significantly reducing the CPU time needed to solve the inverse problem. The presented results show the significant reduction in the CPU time with the proposed method when compared to the standard method.
The rest of this paper is structured as follows. In Section 2 the proposed algorithm is presented and the implementation to a finite and infinite number of degrees of freedom is explained. To verify the accuracy of the proposed algorithm, numerical simulations are performed in Section 3. Furthermore, lab experiments are also used to study the performance of the method which helps to establish the reliability of the proposed algorithm. To this end in Section 4 a prototype software is developed in LabVIEW (National Instruments Ltd.) and used for the lab experiments in two tests where experimental and numerical results are compared to illustrate the efficiency of the proposed algorithm. We finish the paper with some concluding remarks.

Problem Formulation
For a given n-DOFs system subjected to a dynamic excitation u(t) such as shown in Figure 1, applying Newton's second law results into the following system of n-equations written in matrix form where p(t),ṗ(t) andp(t) are the displacement, the velocity and the acceleration vectors, respectively. The matrices M, C and K are the mass, damping and stiffness matrices, respectively. The matrix B u is composed of zeros and ones and is used to apply the load u(t) only at the excitation points. In this study damping is assumed to be proportional to the stiffness. Furthermore, we assume the system is linear and time-invariant while the position of the forces to be identified is determined. It should be noted that the method we propose in this work may also be applied to other discrete systems that can be written in a fashion similar to Equation (1). The state-space equation of 2n-dimensions can be obtained by rewriting the equation of motion (1) such as where A c is the system matrix while B c is the control matrix and are is defined as We divide the time domain into a set of sampling points so that t = t 0 , t 1 , t 2 , . . . , t k , . . . where a uniform time interval δt = t k+1 − t k is used for the divisions. Assuming a reasonably small δt it is possible to approximate the excitation u(t) as a constant within any time interval t k+1 , t k . Thus, the discrete state-space Equation (2) can be rewritten as where x k+1 and x k represent the structural state vector at the time instants (k + 1)δt and kδt, respectively. The constant external excitation in the time interval t k , t k+1 is u k , while A is the system matrix of the discrete state-space equation and B is the control matrix of the discrete state-space equation. The two matrices are given by A = e A c δt (4) We also define the measurement vector y k+1 at the time instant (k + 1)δt as The entries of the matrices H and D are related to the taken measurement which can be displacement, velocity or acceleration. In this paper we use accelerometers to evaluate the acceleration data. Hence, we write the two matrices H and D as where H 0 is a position matrix composed of 0 and 1 with 1 corresponding to the degrees of freedom where a measurement is taken. In this paper we propose an algorithm to identify the excitation load u(t) at each interval in real time as the measurement y k+1 is collected.

Online Force Identification Method
The key point of the proposed method is to build a recursive algorithm for the dynamic load identification in the time domain. The system state is described as a series of time steps that are updated based on the measurement taken and the state unbiased estimation. We then identify the dynamic load using the weighted least square method to minimize covariance unbiased estimation. The state update, the load identification and the measurement update are achieved using the following three equations: The matrices J k and K k are needed for unbiased estimates of the system state i.e., x k|k and the load estimatesû k . The entries of J k are determined using the weighted least square method when identifying the load while those of K k are evaluated within the process of the state unbiased estimation of the minimum variance.
Based on Equations (9) and (3), the error and covariance matrix of x k|k−1 can be written as The covariance matrix of response positions P x k−1|k−1 and P ux k−1 and of excitation positions P u k−1 and P xu k−1 are given by We define the residual vector as y k = y k − Hx k|k−1 (14) where the measurement y k including some error, is defined by By substituting Equation (15) back into (14) we get Substituting Equation (14) into (10)û where x k|k−1 is unbiased estimation, E(v k ) = 0 and E( y k ) = DE(u k ). In addition,û k is unbiased estimation then J k D = I. We also define R k as Using Cholesky decomposition we can rewrite R k as Multiplying Equation (16) with S −1 k we get Thus, the minimum-variance approximation of u k based on the least square method is and the matrix J k is then given by We substitute Equation (10) into (11) to get the state vector x k|k where for convenience we write L k = K k (I − DJ k ). From Equations (23) and (16), we obtain Since x k|k is an unbiased estimation we can assume L k D = 0 and the variance P x k|k of x k|k can be evaluated as Furthermore, because x k|k is unbiased estimation of the minimum variance we aim to minimize P x k|k or its trace tr(P x k|k ), hence, we differentiate tr(P x k|k ) and then search for the solution when the derivative is equal to zero so that while Substituting Equation (27) into (25) we have From Equation (24) and using L k D = 0 we obtain which can now be solved. The previous steps are summarized in Algorithm 1.
2: Dynamic load identification 3: We assume here that the structure parameters including the mass, damping and stiffness matrices are given. The schematic diagram in Figure 2 shows the input and the output of the proposed method. It should be stressed that the response data is updated at the time intervals δt, hence, the dynamic load is identified at the same intervals.

Infinite Number of Degrees of Freedom and Modal Truncation
The computations involved in the above algorithm increases exponentially with the increase in the number of degrees of freedom. Hence, depending on the computational power the above algorithm may not achieve real-time. Here it should be noted that by real time we mean that the CPU time needed to evaluate the dynamic load is an order of magnitude shorter than the time interval δt. Although the algorithm can easily achieve that for relatively small number of degrees of freedom, but this becomes more challenging when studying a complex structure. In this subsection we aim to introduce some changes so that the algorithm will remain efficient even at a large degree of freedom. To this end the modal analysis theory is used to achieve the modal transformation. Hence, the n-order dominant modal characteristic of a structure is used to replace the full modal characteristic. Thus, the computations can be significantly reduced and higher efficiency can be achieved so that the algorithm can still identify the dynamic load in real time even for complex structures. Figure 3, shows a simply supported beam subjected to dynamic load at a given position l a . The beam length is l and the dynamic load functions is u(t). Using the modal transform method, the response of the beam can be given as where q(t) is the modal displacement vector and w = [ϕ 1 , ϕ 2 · · · ϕ n ] the modal matrix. Substituting Equation (31) into (1), we get Since the mode shapes are orthogonal we can write where l = diag(λ 1 , λ 2 · · · λ n ) and λ i = ω 2 i is the ith-order eigenvalue. Then Equation (32) can be simplified asq (t) + gq(t) + lq(t) = w B u u(t) (35) where g = w Cw is modal damping matrix where we assume proportional damping in the system i.e., g = diag(2ξ 1 ω 1 , 2ξ 2 ω 2 , · · · 2ξ n ω n ) with ξ i being the damping ratio of the ith-mode shape.
After transforming the continuous system to a discrete one the approach proposed in the previous section can now be applied where in the modal space the state vector can be defined as while the system A c and the control B c matrices are rewritten as Again, if the measured response is acceleration we can obtain the matrices H and D by Once the matrices A c , B c , H and D are obtained for the continuous system the online-identification algorithm of dynamic loads can be used in the same way as before. But again it should be noted that the number of critical mode shapes to be considered will affect the speed of computations in achieving real time. Basic criteria are often considered for the truncation of the number of mode shapes in vibrations [48]. The natural mode shapes with frequencies close to the excitation frequencies must be considered. In general the first few natural frequencies dominate the response of a structure. Obviously, considering more mode shapes increase the computation accuracy. However, adding more frequencies will also increase the computational cost. Therefore, it is important to identify the number of the critical mode shapes that must be considered to achieve a given accuracy before implementing the algorithm. This will depend on the complexity of the structural as well as the location and the frequencies of the applied excitation. In general for a simply supported beam similar to the one in Figure 3 the first four to eight mode shapes are enough to represent its dynamic behavior.

Numerical Experiments
In this section, we aim to evaluate the efficiency and the accuracy of the proposed algorithm using a set of numerical experiments. Two methods are used to evaluate the results accuracy, namely the peak relative error method (PREM) and the angle cosine method (ACM). For the PREM we use where X(i) represents the exact value while Y(i) is the identified value. The method is only concerned with the error at the maximum value. However, the phase shift is not captured in this method. The error in the ACM is defined as Unlike the previous error measure this method helps to evaluate the global error. The ACM is often used for measuring the error in a waveform signals such as the dynamic loads considered in this paper. The signal similarity is captured by correlation calculation of one-dimensional vectors. The upper limit for this measurement is one. The closer the value of s(X, Y) to one the better the identification results are.
For the numerical experiments we consider a simply supported beam. The beam length is l = 1 m and the cross section is w = 0.05 m×d = 0.005 m 2 . For the material properties we consider steel: the elasticity modulus E = 2.06 × 10 11 Pa, the Poisson's ratio ν = 0.3 and the density ρ = 7900 kg/m 3 . The beam is meshed into 50 8-noded uniform hexahedron elements. The excitation force F(t) is applied at the bottom node between elements 21 and 22. The element numbers and the force location are indicated in Figure 4. The finite element software Patran [49] is used to obtain the beam response under different type of excitations always applied at this node. Using the response measured at the excitation point we then use the proposed algorithm to evaluate the dynamic load. Figure 4 shows the considered finite element. In all the following computations the time interval is δt = 0.001 s and the mode truncation order is six. The choice of this sampling interval meets the Nyquist-Shannon sampling theorem where 0.001 is less than half of the period 0.0024.

Multiple Impacts at Fixed Intervals
First, we consider an impact load applied as half sine wave with a duration of 0.02 s and an amplitude 1N. Four impacts are applied at 3 s intervals. The application time instants are t 1 = 0.1 s, t 2 = 3.1 s, t 3 = 6.1 s and t 4 = 9.1 s as shown in Figure 5. The measured response is obtained using the finite element model. The response is fed into the algorithm without applying any noise and the accurate mode shapes are considered. Next, the 5% zero-mean Gaussian noise is added to the measured response as well as a 5% deviation is added to all the modal parameters so that the used mode shapes are inaccurate. The four possible combination of the clean and contaminated data are considered for the algorithm input. The contaminated data is used to study the stability of the algorithm. Figure 6a shows the identified results using the proposed algorithm. The figure shows that the peak load value is captured with good accuracy (relative error is 2.91%),which is the most critical value to be identified in the case of impact load. However, the figure also shows some noise in between the peaks which has resulted in inflating the errors to some extent. The PREM and ACM errors for all these cases are listed in Table 1.   The results suggest that the proposed online-identification method can efficiently reconstruct the dynamic loads. Due to the limited number of mode shapes we can observe some error in identifying the exact dynamic load. If we look into the results in Table 1 we can see that both the PREM and the ACM errors are reasonable. Also, the identified load pattern is well in agreement with the exact load pattern as can be seen in Figure 6. Clearly, when we use polluted data the error increases but this increase in the error remains in line with the inaccuracy in the data i.e., introducing 5% noise into the input leads to a proportional increase in the identification error. This suggests that the algorithm is stable and can produce meaningful results even when the data is contaminated with noise.

Multi-Frequency Sinusoidal Load
Next we aim to test the proposed algorithm for identifying a periodic load. We redefine the excitation force as F(t) = 0.5 sin(40πt) + 1.5 sin(60πt). All the remaining parameters of the model are kept unchanged. However, we measure the beam response at node A which is indicated in Figure 4. Unlike the previous case here we collect the beam response in a point that is different from the excitation point. We use the proposed algorithm to reconstruct the dynamic load in time steps. Again we consider clean and polluted data as an input for the four different cases considered in the previous test.  Table 2 the corresponding errors in the four considered cases. As can be expected the best approximation is achieved in the case where clean data is used i.e., Figure 7a where the PREM and the ACM errors are about 0.38% and 0.9741, respectively. This relatively small error was achieved using only six mode shapes in the approximation. The 5% Gaussian white noise introduced into the measured response has a small impact on both types of errors. However, a 5% error in the modal properties cause an order of magnitude increase in the PREM error or even two orders of magnitude when added to the Gaussian white noise. But it should be also noted that despite the significant increase in the error when identifying the peak load but the global error i.e., ACM remains reasonably small and the identified wave pattern matches with high accuracy the actual sinusoidal load.

Experimental Validation
In this section, we aim to validate the proposed algorithm using physical experiments. We run the same experiments as before but using a physical beam under dynamic loading. To this end it is important to build a signal acquisition system and integrate it with the algorithm. In other words we need to build both experimental hardware and software to validate the proposed algorithm in a physical experiment. The hardware include the sensor and its connections and a signal acquisition module, a vibration exciter and a computer as well as the beam. The software is implemented in LabVIEW (National Instruments Ltd.) where beside the proposed algorithm we also use the modules: response acquisition and display, modal data reading, online identification of dynamic load and data storage and recall. The components of the experimental work are displayed in Figure 8. The algorithm implementation on LabVIEW (National Instruments Ltd.) is shown in Figure 9 and the source code is given in the Appendix A. After starting the system, the parameters are first initialized and then the software waits for operation command. Once the system is operational it initiates the mode data reading and Kalman Filter. Next, the data collection and the load identification begins. The response data collected, and real-time identified load is then displayed step-by-step on the user interface. To realize the software real-time load identification function, it is necessary to perform four different tasks simultaneously: uninterrupted response data acquisition, Kalman Filter iterative identification of the load, interface display and continuous force excitation. This puts a significant computational burden on the computing hardware. To deal with this it is important to develop a parallel implementation for the system. The programming language of LabVIEW, is inherently multi-threaded which helps to design a suitable implementation.
We also show in Figure 10 the load identification interface which shows the needed parameters including mode truncation order, model node number, observation node number and excitation point position. Other parameters such as the natural frequencies and damping can be loaded from a saved file. Also, here we validate the implementation using the same two previous cases, namely multiple impacts at fixed intervals and constant frequency sinusoidal load. Figure 11 shows a schematic diagram for both experiments. In Figure 12 we show the experimental settings. The used exciter is JZT-2 while as an excitation source either the power amplifier (HEAS-50) or the force hammer (PCB 086C03) are used. The acceleration sensor (PCB 356A33) and acquisition board (NI USB 4431) are used to collect the response requested by the online dynamic load identification software developed to identify the load acting on the beam. It should be noted that the weight of the acceleration sensor is negligible relative to the weight of the beam. The experimental modal analysis is curried out using a force hammer to excite the beam at several points. Then the modal information such as natural frequency and modal damping are obtained automatically by the modal analyzer (M + P Smart Office).

Multiple Impacts at Fixed Intervals
Now we aim to validate the proposed algorithm implementation using an experimental approach. Again, we consider a simply supported beam of the lab. In this work the beam length, width and thickness are a = 0.695 m, b = 0.04 m and h = 0.008 m, respectively. The model beam has a rectangular cross section w = b × h = 0.00032 m 2 . The material properties are: density ρ = 7900 kg/m 3 , modulus of elasticity E = 2.06 × 10 11 Pa and Poisson's ratio ν = 0.30. The time interval is δt = 0.001 s. The sampling rate is 1024 Hz. It should note that the length and the width of the considered beam has no impact on the algorithm performance i.e., we expect the same conclusions when repeating the experiment with beams of other dimensions. The first seven mode shapes and their properties are evaluated experimentally and listed in Table 3. To verify the implementation, we impose an impact load and try to identify it in real time using the data collected by the accelerator sensor. The load identification parameters in the software are shown in Figure 13. The actual load is then compared to the identified load in real time. The load identification results of the first and second impacts are plotted in Figure 13. The plots zoom on the results within 0.1 s of the impact. The algorithm recovers with good accuracy the peak value of the first impact. In general the identified load matches the actual load up to the peak value of the impact.
Thereafter the identification results start to fluctuate due to the observation noise and model error before it is stable again within 0.02 s. Similar observation can also be made for the second impact also plotted in the figure. Again, the peak value is recovered with good accuracy and some fluctuations are observed beyond the peak value. The errors with the PREM and the ACM are, respectively, 4.15% and 0.5839. The relatively large global error can be attributed to the fluctuations in the recovered load. Furthermore, the error in the measured damping compared to the actual damping can also contribute to these oscillations.

Constant Frequency Sinusoidal Load
In the last experiment we aim to validate the proposed implementation when recovering a multi-frequency sinusoidal load. This loading case will activate fewer mode shapes compared to the previous case. The same beam is considered again, and all the parameters are kept unchanged from the last experiment with only replacing the load. The input parameters are shown in Figure 14. Also in the same figure we plot the identified load compared to the actual load. The plots also show the acceleration recovered by the sensor and used as an input for the algorithm. In the initial stage we can observe a small drift in the amplitude and the phase between the identified and the actual. This is stabilized in around 0.5s and soon after two cycles of loading. The errors obtained in this case are 2.64% and 0.9381 with the PREM and the ACM, respectively. This reflects the accuracy in identifying the load.

Comparison with Newmark Method
After validating the approach using both numerical and experimental tests we now want to evaluate the approach performance relative to a standard inverse method. We choose Newmark method which is often considered to be an industrial standard method that is used to evaluate the dynamic loads offline. For the comparison purpose we consider a simply supported beam with l = 1 m length and w = 0.05 m×d = 0.01 m 2 cross section and the material properties: elasticity modulus E = 2.1 × 10 11 Pa, Poisson's ratio ν = 0.3 and density ρ = 7800 Kg/m 3 . The beam is meshed using 10 uniform 8-noded hexahedron elements. The applied excitation is F(t) = 20 sin(80πt) which is located at the bottom node between elements 3 and 4 while the beam response is measured at node A all of which are indicated in Figure 15. The remaining parameters are kept unchanged from previous tests. The computations were run on Windows PC with an Intel Core(TM)i7 3.6 GHz CPU and 32 GB of RAM. It should be noted that a similar example is also published in [50]. For the comparisons we again consider the input data without any noise. We then add a Gaussian noise with zero-mean and signal-to-noise ratio of 5%. We consider a third case where a 5% deviation is added to all the modal parameters and a fourth case where both the noise and the deviation are considered together. The contaminated data will help to compare the stability of the proposed algorithm compared to Newmark method. Figure 16a-d show the results of the proposed algorithm plotted against the exact dynamic load as well as those obtained with Newmark method. The PREM and ACM errors for all these cases are listed in Table 4. The results show that both methods can recover the dynamic load with good accuracy when using clean data. Furthermore, the proposed algorithm has an improved accuracy compared to Newmark method. This is true for the global error PREM as well as the error in the peak value ACM. However, the proposed algorithm was able to evaluate each time step within 0.4 ms compared to 8ms with Newmark. This is equivalent to a 95% reduction in the CPU time with the proposed algorithm. Moreover, since the considered time step is 1ms it is possible to evaluate each time step in real time as soon as the input data is measured whereas the CPU time with Newmark method is eight times longer than the time step. The results in Figure 16b shows that the accuracy of Newmark method is significantly affected by adding the Gaussian noise where the results become meaningless as can be seen in both PREM and ACM errors. However, the proposed method still produces meaningful results and despite the relatively high PREM error, the peak value is still measured with good accuracy as reflected by the ACM error. The Newmark method seems to be less sensitive when using inaccurate modal properties, but the proposed method remains to be more accurate. Again when considering a combination of polluted data with inaccurate modal properties the proposed method still produce meaningful results which is reflected in the errors as well as Figure 16d while Newmark method again fail to produce useful results. This results suggest that the proposed algorithm is more precise and more stable against noise data when compared with Newmark method. This improved performance is achieved with a major reduction in the CPU.

Conclusions
In this paper, we propose an algorithm to identify in real time the dynamic forces acting on a vibrating structure. The algorithm that is based on Kalman filter, consists of three steps. First, the system and the control matrices are built for the vibrating structure using the vibration differential equation for a discrete system or the mode truncation method for a continuous system. Next, the weighted least square method is used to minimize covariance unbiased estimation. Then a state vector is updated based on acceleration measurements. By minimizing the prior estimate error covariance, the excitation force is estimated. This is repeated at a fixed time interval where the applied load is updated in the time domain. As an application the proposed algorithm is used to identify the dynamic load applied on a simply supported beam. The performed numerical simulations show that the proposed method can be used to find accurate estimate for the applied force even when the measured data or the model parameters are not very accurate. To validate the proposed algorithm experimentally, a LabVIEW (National Instruments Ltd.) implementation is presented and discussed. The proposed algorithm is verified through two different types of dynamic loads. The experimental results show a good potential for the proposed approach in identifying the dynamic load applied on a structure in real time. Moreover, when compared to a standard method the proposed algorithm is more accurate and significantly more stable. The results show 95% reduction in the CPU time required with the standard method. The work presents one of the first attempts in real time dynamic load identification the and can have several future applications.

Conflicts of Interest:
The authors declare no conflict of interest.