Input Forces Estimation for Nonlinear Systems by Applying a Square-Root Cubature Kalman Filter

This work presents a novel inverse algorithm to estimate time-varying input forces in nonlinear beam systems. With the system parameters determined, the input forces can be estimated in real-time from dynamic responses, which can be used for structural health monitoring. In the process of input forces estimation, the Runge-Kutta fourth-order algorithm was employed to discretize the state equations; a square-root cubature Kalman filter (SRCKF) was employed to suppress white noise; the residual innovation sequences, a priori state estimate, gain matrix, and innovation covariance generated by SRCKF were employed to estimate the magnitude and location of input forces by using a nonlinear estimator. The nonlinear estimator was based on the least squares method. Numerical simulations of a large deflection beam and an experiment of a linear beam constrained by a nonlinear spring were employed. The results demonstrated accuracy of the nonlinear algorithm.


Introduction
Advanced structural health monitoring is generally regarded as a vital technology for the next generation of aeronautical and space systems [1]. This technology is aimed at preventing catastrophic structural failures and is comprised of three facets: (a) determination of stresses and deformations of structural components; (b) estimating input forces; and (c) detection of critical damage mechanisms such as cracking, delamination, and corrosion. Estimating dynamic input forces of engineering structures accurately has great significance for structural health monitoring, fatigue analysis and life estimation, and is regarded as the premise of structural design and optimization. Damage tolerance considerations determine the design of composite components, and forces occurring in service conditions can induce non-visible damages in structures. Besides, additional uncertainties may arise, as structural architectures behave in complex ways under the action of thermo-mechanical forces. By knowing forces distribution, we can adjust structure layout and enhance material performance to ensure safety. Moreover, the presence of damage can modify load paths in difficulty predictable ways. By knowing force distribution, we can also assess damage as there is usually structural damage in the location of the forces mutation. In engineering applications, input forces are quite hard or even impossible to be measured directly. Some examples are the measurements of airplane wing deflection, blade shape changes of windmill or helicopters, tool-tip displacement of line boring machines, etc. In these cases, it is helpful to use attached sensors to measure the responses of structures. Estimating dynamic input forces by using responses of structures is an inverse problem.
In mathematics, inverse problems are typically ill-posed and are difficult to solve. For inverse problems, a little measurement errors can cause large estimation errors. So input forces estimation method is the key to improve accuracy.
In summary, input forces estimation algorithms [2][3][4][5][6][7] can be divided into time domain algorithms and frequency domain algorithms. Time domain algorithms are much more valuable than frequency domain algorithms, as time domain algorithms can accomplish estimation in real-time. In dealing with nonlinear beam systems, there is real-time change in the stiffness matrix, which means that frequency domain algorithms cannot be used at this situation. In the process of the input forces estimation, suppressing noise is the key to improve accuracy, as small measurement errors can cause large estimation errors in inverse problems. For most practical engineering applications, noises can be simplified as white Gaussian noise. This means that the method of suppressing white Gaussian noise need to be studied first. Ma et al. [8,9] proposed a method to estimate input forces of linear beam systems by combining Kalman filter with a recursive method. In his work, Kalman filter was used to suppress noise, and residual innovation sequences, a priori state estimate, and innovation covariance generated by Kalman filter were used to estimate input forces by using a least-squares method. However, these studies were associated with linear beam systems, and estimating dynamic input forces of nonlinear beam systems have not been studied. As nonlinear beam systems are widely used in engineering applications, this work focus on estimating input forces for nonlinear beam systems.
In the work, with the system parameters determined, the magnitude and location of input forces in nonlinear systems can be estimated by using SRCKF and a nonlinear estimator. The nonlinear estimator is based on least squares method. According to the second order dynamic system and measuring principle, the state equations and measurement equations of the state-space model are established. The Runge-Kutta fourth-order algorithm is employed to discrete the state equations. SRCKF is used to suppress noise, and the residual innovation sequences, a priori state estimate, gain matrix and innovation covariance generated by SRCKF are employed to estimate the magnitude and location of input forces by using a nonlinear estimator. To verify the effectiveness of this estimation method, numerical simulations of a large deflection beam and experiment of a linear beam constrained by a nonlinear spring are employed.

Problem Formulation
There are three steps to estimate the location and magnitude of input forces. First, the state equation and measurement equation of the state-space model are discretized by using a Runge-Kutta method. Second, SRCKF is used to suppress white noise. Finally, residual innovation sequences, a priori state estimate, gain matrix and innovation covariance generated by SRCKF are used to estimate the location and magnitude of input forces.

Discretization of a Nonlinear System
In the paper, the dynamic parameters of nonlinear beam structure were known, and we could then construct the state-space model. With unknown parameters of structure in engineering applications, we could estimate dynamic and static parameters of state-space models according to Refs. [10][11][12][13]. The Runge-Kutta methods are widely used to discrete continuous-time system, and the most well-known member Runge-Kutta fourth-order algorithm (RK4) has the advantages of high precision, convergence and stability. Therefore, RK4 is employed by the paper. The nonlinear, continuous-time model can be described by: . .
where state vector X(t) = [position; velocity] and observation vector Z(t) represent measured dynamic responses; vector F(t) = [0; F]; f 0 (·) and h 0 (·) are nonlinear functions with respect to X, F and t; F is the forces vector; w(t) and v(t) represent continuous-time white noise processes. For a step-size ∆T > 0 and an initial value X k−1 , the state Equation (1) can be discrete as follows: for k = 1, 2, 3, 4, . . . , using: Here, X k is the state approximation, and the value X k is determined by the value of X k−1 plus the weighted average of the increments (b 1 , b 2 , b 3 , and b 4 ), where each increment is the product of the sample interval.
The discrete model is described by: where vector w k represents the process white noise, Q represents covariance matrix and δ kl is the Kronecker deltas.
where vector v k represents the measurement white noise, R represents the noise covariance matrix, R v = σ 2 , σ is the standard deviation of the measurement noise. The vectors w k and v k are mutually uncorrelated.
The initial value S 0|0 of the square-root factor of the error covariance matrix can be computed by the Cholesky decomposition, where chol[·] is the Cholesky factorization.

Time Update
(1) Calculate the cubature points: wherex k/k is the prior estimated state. ξ i = m 2 [1] i , [1] i is the ith column of the matrix [I (−1)I].
(2) Calculate the propagated cubature points: (3) Calculate the predicted state:x (4) Calculate the square-root factor of prediction error covariance: where S Q,k is obtained by Cholesky factorization for Q k = S Q,k S T Q,K . tria(·) is a matrix triangularization algorithm which can generate a lower triangular matrix:

Measurement Update
(1) Calculate the cubature points: (2) Calculate the propagated cubature points: (3) Calculate the predicted state:ẑ (4) Calculate the square-root of the innovation covariance matrix: is obtained by Cholesky factorization for R k+1 = S R,k+1 S T R,k+1 and R k+1 is the noise covariance matrix. tria(·) is a matrix triangularization algorithm which can generate a lower triangular matrix.

The Nonlinear Estimator
By applying residual innovation sequences, a priori state estimate, gain matrix and innovation covariance generated by SRCKF, input forces can be estimated by using a nonlinear estimator from the response values (displacement, velocity, or acceleration). The inverse estimation method consists of two parts: SRCKF with no input forces terms, and a nonlinear estimator. In the nonlinear estimator, the first-order Taylor series expansion is used to arrive at the estimated state valuex k/k−1 , and a least squares method is used to estimate forces. The detailed derivation of the nonlinear estimator can be found in Appendix A. The simple equations of the nonlinear estimator are as follows: where f (·) and h(·) represent nonlinear functions of the discrete system, P zz,k/k−1 represents the innovation covariance matrix, K k represents the gain matrix, B s (k) and M s (k) represent the sensitivity matrices, Z k represents the innovation matrix, K b (k) represents the correction gain for updatinĝ F(k), P b (k) represents the error covariance,F(k) represents the estimated input vector, and γ is a fading factor. The procedures for the nonlinear method are summarized as Figure 1.

Simulation Model
In structural dynamic analysis, the slender beam may exhibit geometrically nonlinear behaviors when it undergo large deformation. In the paper, a large deformation beam is picked as the model, and the equation of motion can be described as follows:

Simulation Model
In structural dynamic analysis, the slender beam may exhibit geometrically nonlinear behaviors when it undergo large deformation. In the paper, a large deformation beam is picked as the model, and the equation of motion can be described as follows: where M is the mass matrix, C the damping matrix, K the stiffness matrix, X the displacement vector, and F the equivalent nodal force vector.
It is assumed that the beam element is two-dimensional and each node has three degrees of freedom (two translational and one rotational). Beam element is shown as Figure 2.

Simulation Model
In structural dynamic analysis, the slender beam may exhibit geometrically nonlinear behaviors when it undergo large deformation. In the paper, a large deformation beam is picked as the model, and the equation of motion can be described as follows: where M is the mass matrix, C the damping matrix, K the stiffness matrix, X the displacement vector, and F the equivalent nodal force vector.
It is assumed that the beam element is two-dimensional and each node has three degrees of freedom (two translational and one rotational). Beam element is shown as Figure 2. According to Refs. [16][17][18], the element linear stiffness matrix , nonlinear stiffness matrix and consistent mass matrix can be expressed in the form: According to Refs. [16][17][18], the element linear stiffness matrix K e L , nonlinear stiffness matrix K e N and consistent mass matrix M e can be expressed in the form: where ρ is mass density, A is the area of cross-section, L is the beam element length, E is Young's modulus of elasticity, I is the moment of inertia of the cross-section, and u 1 u 4 are the axial deformation of the first node and second node. In converting the second order dynamic system to the state-space model, the state equation and measurement equation can be written as: . where: is a nonlinear function with respect to Y. H is a measurement matrix and Z(t) represents the measurement values vector.
Equations (29) and (30) are discretized using RK4, and the discrete model can be described by: where vector w k represents the process white noise, Q represents the covariance matrix, and δ kl is the Kronecker deltas.
where vector v k represents the measurement's white noise, R represents the noise covariance matrix, R v = σ 2 , and σ is the standard deviation of the measurement noise. The vectors w k and v k are mutually uncorrelated.
Considering a five-element beam, the parameters of the beam are: Elastic modulus E = 7.2 × 10 10 N/m 2 ; density ρ = 2.7 × 10 3 kg/m 3 ; beam length l = 1 m; cross section S = 0.1 m × 0.01 m; sampling interval ∆T = 0.001 s. The system responses (displacements and rotations) are obtained by RK4, and the responses with white noise are employed as the measured dynamic responses. Thus, the magnitude and location of input forces can be estimated in real-time from dynamic responses (displacements and rotations).
The initial parameters of the estimation system are generally listed as follows: x 0 = 0 30×1 , P 1 = I 30×30 , P 2 = 0 30×30 , M s = 200 × I 30×30 , P b = 200 × I 30×30 , γ = 0.69. To verify the effectiveness of this estimation method, the root-mean-square error (RMSE) method is used to measure the errors between the estimated forcesF i and the exact forces F i . The RMSE is computed as follows:

Simulation Results
(1) For the large deflection beam model, a sinusoidal input force, a rectangular input force, and a triangular input force are estimated, respectively. Figures 3-8 plot the results of magnitude estimation. For the method proposed in the paper, the forces of total degrees of freedom can be estimated, so we can estimate the location of force if single force is applied. In the simulation, force is applied at the sixth degree of freedom with a total of 15 degrees of freedom. The results of estimation of 15 forces at 15 degrees of freedom are plotted at Figures 9-11.
(2) The performance (the mean errors and RMSE values) of estimation method with Q w at 1 × 10 −4 , 1 × 10 −6 and σ at a constant 1 × 10 −8 are presented in Table 1. estimated, so we can estimate the location of force if single force is applied. In the simulation, force is applied at the sixth degree of freedom with a total of 15 degrees of freedom. The results of estimation of 15 forces at 15 degrees of freedom are plotted at Figures 9-11.
(2) The performance (the mean errors and RMSE values) of estimation method with at 1 × 10 −4 , 1 × 10 −6 and at a constant 1 × 10 −8 are presented in Table 1.   estimation. For the method proposed in the paper, the forces of total degrees of freedom can be estimated, so we can estimate the location of force if single force is applied. In the simulation, force is applied at the sixth degree of freedom with a total of 15 degrees of freedom. The results of estimation of 15 forces at 15 degrees of freedom are plotted at Figures 9-11.

Discussions of the Simulations
(1) The SRCKF is famous for strong stability and higher precision, and only needs recent measurement data and the previous estimated value to estimate input forces, and so the proposed method could save computer memory, reduce computational burdens, and improve system robustness.
(2) Figures 5-8 show that the estimated input forces rapidly converge to exact input forces with non-zero initial state, but with a large initial estimation errors. Table 1 show that the estimation performance of the sinusoidal input force are better than that of the rectangular input force and triangular input force, and this is because abrupt changes of input force will cause large errors. Figures 4, 6 and 8 show that the system have good stability with large system noise and measurement noise. From Figures 9-11, we can conclude that the location estimation of forces have a good performance. For three types of forces, the estimation performance of the sinusoidal input force and the triangular input force are better than that of the rectangular input force.
(3) The results in Table 1 show that the proposed estimation method has good capabilities to suppress noise. The mean errors for sinusoidal input force, rectangular input force and triangular input force are close to zero. The RMSE values for sinusoidal input force, rectangular input force, and triangular input force are less than 0.46%, 6.61%, and 5.27%, respectively.

Experimental Model and Measurement Principle
Considering a cantilever with a nonlinear spring stalled at the end node, the finite element model and the FBG (Fiber Bragg Grating) sensing network are shown in Figure 12. In the experiment, the strain values of FBG sensor network were used as observed values. The parameters of the beam are: Elastic modulus E = 6.89 × 10 10 N/m 2 ; Density ρ = 2.69 × 10 3 kg/m 3 ; Beam length l = 0.48 m; Cross section S = 0.03 m × 0.003 m. The performance of the nonlinear spring is plotted in Figure 13. There were six measuring points of FBG on the beam along its center line which were employed to record the beam's surface strain simultaneously. The distance between two consecutive sensors was about 9.15 cm. Figures 4, 6, and 8 show that the system have good stability with large system noise and measurement noise. From Figures 9-11, we can conclude that the location estimation of forces have a good performance. For three types of forces, the estimation performance of the sinusoidal input force and the triangular input force are better than that of the rectangular input force.
(3) The results in Table 1 show that the proposed estimation method has good capabilities to suppress noise. The mean errors for sinusoidal input force, rectangular input force and triangular input force are close to zero. The RMSE values for sinusoidal input force, rectangular input force, and triangular input force are less than 0.46%, 6.61%, and 5.27%, respectively.

Experimental Model and Measurement Principle
Considering a cantilever with a nonlinear spring stalled at the end node, the finite element model and the FBG (Fiber Bragg Grating) sensing network are shown in Figure 12. In the experiment, the strain values of FBG sensor network were used as observed values. The parameters of the beam are: Elastic modulus = 6.89 × 10 (N/m ); Density = 2.69 × 10 (kg/m ) ; Beam length = 0.48 m; Cross section = 0.03 m × 0.003 m. The performance of the nonlinear spring is plotted in Figure 13. There were six measuring points of FBG on the beam along its center line which were employed to record the beam's surface strain simultaneously. The distance between two consecutive sensors was about 9.15 cm.   The relationship between wavelength shift of FBG sensors and strain values is showed in Equation (33). Supposing with no temperature change, the strain value can be computed by measuring the wavelength change according to Equation (34): where ∆ is the wavelength shift, is the axial strain, is the thermal expansion coefficient, is the thermo-optical coefficient, ∆ is the temperature change, and is the effective photo-elastic The relationship between wavelength shift of FBG sensors and strain values is showed in Equation (33). Supposing with no temperature change, the strain value can be computed by measuring the wavelength change according to Equation (34): where ∆λ is the wavelength shift, ε is the axial strain, α is the thermal expansion coefficient, ξ is the thermo-optical coefficient, ∆T is the temperature change, and P e is the effective photo-elastic coefficient [19][20][21]. For a Euler-Bernoulli beam with n-element, the relationship between strain values and nodal degrees of freedom can be described as follows: is the shape function; l is the length of the beam element; ξ = x/l, x is the location of the FBG in element; h is the thickness of the beam; w is nodal displacement; θ is nodal rotation. [ε 1 , ε 2 , · · · ε 2n ] T is strain values vector. For the experiment with six FBG sensors pasted on beam, observation matrix can be describes as follow:

Experimental Procedure and Results
The layout of experiment is plotted in Figure 14. FBG (Fiber Bragg Grating) interrogation system (SM130 (SonMicro Elektronik, Mersin, Turkey)) is used for measuring the dynamic strains, and an electrodynamics shaker is employed for the excitation. Input force is applied at the beam end, and the magnitude and location of input force are estimated by FBG sensor network. In the experiment, a force sensor is stalled between the exciter and beam end, and the measured force is used as an exact value to verify the practicability of the proposed method. In the process, a NI cDAQ-9174 module (National Instruments Corporation, Austin, TX, USA) and LABVIEW software (National Instruments Corporation, Austin, TX, USA) were used to acquire the signal. The sampling frequency was set as 100 Hz, and the experimental time is was two seconds. The results of input force estimation were plotted in Figure 15. a force sensor is stalled between the exciter and beam end, and the measured force is used as an exact value to verify the practicability of the proposed method. In the process, a NI cDAQ-9174 module (National Instruments Corporation, Austin, TX, USA) and LABVIEW software (National Instruments Corporation, Austin, TX, USA) were used to acquire the signal. The sampling frequency was set as 100 Hz, and the experimental time is was two seconds. The results of input force estimation were plotted in Figure 15.  a force sensor is stalled between the exciter and beam end, and the measured force is used as an exact value to verify the practicability of the proposed method. In the process, a NI cDAQ-9174 module (National Instruments Corporation, Austin, TX, USA) and LABVIEW software (National Instruments Corporation, Austin, TX, USA) were used to acquire the signal. The sampling frequency was set as 100 Hz, and the experimental time is was two seconds. The results of input force estimation were plotted in Figure 15.

Discussions
(1) Experimental results showed that the estimated input forces had a little amplitude error. The amplitude error was mainly produced by insufficient number of sensors. In engineering applications, distributed optical fiber sensors network can solve the deficiency of sensors installation, where the distance between two consecutive sensors can be less than 1 cm.
(2) The proposed estimation method required the statistical characteristics of noise to be known, as well as an accurate system model. In practice, sometimes the statistical characteristics of noise are unknown, and the system model is inaccurate. In addition, the nonlinear system is easily affected by the model uncertainties in the actual operating environment. For these deficiencies, we employed an adaptive algorithm to estimate the time-varying noise statistics and model uncertainties.

Conclusions
A real-time nonlinear method for estimating input forces is presented in this work. The method used SRCKF to suppress noise and a nonlinear estimator to estimate input forces. Simulations of the large deflection beam system and experiment of a linear beam constrained by a nonlinear spring were applied. Simulation results showed that the mean errors and RMSE values of three types of input forces subjected to the above noise were satisfied. Experimental results showed that the estimated input forces had a little amplitude error, and the estimation method had good stability. Author Contributions: Xuegang Song proposed the methods, performed the experiments and wrote the paper; Dakai Liang gave some suggestions for the improved method; Yuexin Zhang provides important suggestions for improving the paper.

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

Appendix A
Following the nonlinear discrete system as given in: where w(k) and v(k) are white noise sequences, The posteriori state estimate without the exciting force: The posteriori state estimate with the exciting force: where K k is got from SRCKF, and: Define the difference of the two posteriori state estimate as follows: Assume the exciting force begin with the time t n , then: In summary, we get: At the time t n+1 , Equation (A7) becomes: From (A7) we know that ∆X n = 0, so Equation (A8) becomes: Then Equation (A9) becomes: From Equations (A7) and (A10), for k > n, we have: Ignoring Φ k ∆X k−1 , and Equation (A11) becomes: From Equations (A12) and (A11), we have: Assume F k−1 = F k−2 , then: From Equations (A12) and (A13), we have: where: The observed value of the residual sequence with exciting force can be described as: The observed value of the residual sequence without exciting force can be described as: For different values of k, we have: In summary, we get: where: B k = H k (Φ k M k−1 + I)Γ k For k = n + 1, n + 2, . . . , n + l. we have: where: We insert the following term between B T (N + 1) and s −1 (N + 1), and we will get the outcomes: