A novel strategy for force identification of nonlinear structures

Abstract

Dynamic force is the key indicator for monitoring the condition of a mechanical product. These mechanical structures always encompass some nonlinear factors. Most previous studies focused on obtaining the dynamic force of linear structures. Consequently, this study focuses on the nonlinear mechanical structure, and a novel identification strategy is proposed to indirectly identify the excitation force. For the identification strategy, based on a nonlinear state-space model, a force identification equation for the nonlinear structure is built, wherein the transfer matrix consists of coefficient matrices of the nonlinear state-space model, and these coefficient matrices are calculated by a nonlinear subspace identification algorithm. Then, under the generalized cross-validation criterion, the truncated total least squares method is introduced to solve the ill-posed equation to eventually obtain the excitation force of the nonlinear structure. The identification results from two numerical simulation cases and one experimental case illustrate that the proposed identification strategy can stably and accurately identify the excitation force of nonlinear structures.

Introduction

To implement reliability analysis, vibration control, optimal design, and structural health monitoring, it is necessary to obtain accurate dynamic forces acting on the mechanical structure.^1,2 Unfortunately, because of the limitation of installing space and sensing technology, it is always impossible to directly measure dynamic force by the force sensor in most engineering practices. However, compared with measuring dynamic force, the structural dynamic response can be easily measured by a suitable sensor. Therefore, many scholars have devoted time to developing some indirect techniques to indirectly obtain the unknown dynamic force on the basis of the dynamic response and structural characteristics. ³

Indirect force identification is essentially the second kind of dynamic inverse problem, and the key step is to establish a force identification equation and perform the inverse operation of its transfer matrix. Unfortunately, the transfer matrix for the mechanical structure is ill-conditioned. This makes the corresponding force identification problem an ill-posed inverse problem. Consequently, a direct matrix inversion operation cannot be used to solve such ill-posed problems.

In the past several decades, various regularization methods have been developed to solve this problem.^4,5 Qiao and Rahmatalla ⁶ used the generalized cross-validation (GCV) criterion combined with Tikhonov regularization approach to estimate the moving force on an Euler–Bernoulli beam, and they also studied the influence of different regularization matrices on identification results. By identifying the excitation force on a plate, Choi et al. ⁷ studied the influence of regularization parameter selection criterions on the performance of Tikhonov regularization method. This demonstrates that compared with GCV criterion and ordinary cross-validation criterion, the L-curve criterion is more appropriate for obtaining the Tikhonov regularization parameter at a high noise level. Under determining the truncation threshold using the L-curve criterion, Leclere et al. ⁸ developed truncated singular value decomposition (TSVD) method to indirectly identify bearing loads inside an operating engine. Gunawan ⁹ introduced the trust region strategy to enforce the regularization and extended the Levenberg–Marquardt iterative algorithm to solve an unknown impact force. Wang et al. ¹⁰ proposed a novel hybrid conjugate gradient algorithm to iteratively solve multisource force identification problems, and the identification results showed that the stability, robustness, and computational efficiency of the proposed iterative algorithm were superior to those of the traditional Landweber iteration regularization method. Since the profile of the impact force is similar to that of a cubic B-spline, Liu et al. ¹¹ presented a modified regularized cubic B-spline collocation method to identify the impact force. Because of the sparsity of the impact signal, Liu et al. ¹² introduced a fast iterative shrinkage-thresholding algorithm to solve a sparse regularization model developed from the dynamic response and impact force reconstruction equation, so as to synthetically estimate the above two physical quantities. Qiao et al. ¹³ proposed the group sparse regularization approach, which is based on the mixed l_2,1-norm regularization model and an accelerated gradient descent algorithm to identify unknown impact forces. He et al. ¹⁴ combined a non-probabilistic analysis model and TSVD method to identify the load in an acoustic structure interaction system involving non-probabilistic uncertainty. Additionally, He et al. ¹⁵ proposed a novel modified regularization method to identify the random dynamic load. This method is not sensitive to the selection of measuring positions.

In conclusion, various regularization methods have been widely used in engineering practice to successfully solve the problem of indirectly obtaining dynamic loads. However, most of these studies focused on linear structures. In recent years, with the development of nonlinear dynamics including nonlinear numerical algorithms,^16-19 fractal vibration theory,^20-22 and nonlinear system identification,^23-25 nonlinear factors and phenomena in mechanical structures have attracted increasing attention. This means that mechanical structures are always nonlinear because of nonlinearities in material properties and structural components,^26-28 and considering these nonlinear factors can obtain more reasonable and accurate force identification results. These nonlinear factors change the structural characteristics, making it impossible to form the force identification equation as in the previous approach. Furthermore, the nonlinear force caused by the nonlinear factor also reduces the identification accuracy of the excitation force. To identify the non-harmonic periodic excitation force of the nonlinear damped system, Jang et al. ²⁹ introduced Landweber–Fridman regularization and Tikhonov regularization methods to solve an ill-posed Volterra-type nonlinear integral equation in order to obtain a stable identification result. However, this approach requires both the velocity and displacement responses and is only suitable for a single-degree-of-freedom (DOF) system. Ma et al. ³⁰ developed Tikhonov regularization combined with GCV criterion to solve the force identification equation of nonlinear multi-degrees-of-freedom system in the frequency domain. For the frequency-domain method, the identification accuracy near the resonance frequency is low. Additionally, this study only considered errors in the measured dynamic response. However, for nonlinear structures, it is difficult to obtain an absolutely accurate transfer matrix. This means that there are errors in both transfer matrix and response vector.

The indirect force identification problem for a nonlinear structure is studied in this study, and a novel time-domain identification strategy is proposed. The core of the novel strategy is to establish a time-domain force identification equation of nonlinear structure through a nonlinear state-space model, wherein the transfer matrix consists of its coefficient matrices. Therefore, before identifying the excitation force, these coefficient matrices are estimated by nonlinear subspace identification (NSI) algorithm.^31,32 Then, considering that there are errors in both the transfer matrix and dynamic response, the ill-posed force identification equation is solved by the truncated total least square (TTLS) method^33,34 in conjunction with GCV criterion, eventually obtaining the unknown excitation force of the nonlinear structure.

Force Identification Equation Based on Nonlinear State-Space Model

Force identification is the second type of inverse problem in vibration research, and the purpose is to estimate the external load based on structural transfer characteristics and output responses. Therefore, to conduct force identification, the first problem is to determine the relationship between the dynamic response and the excitation force according to structural characteristics, that is, to form a force identification equation. For the linear structure, the convolution of excitation force $f\left(t\right)$ and unit impulse response function $h\left(t\right)$ is equivalent to the dynamic response $y\left(t\right)$ , and it is given by $$y\left(t\right)=h\left(t\right)\otimes f\left(t\right)={\displaystyle {\int }_0^{t}h\left(t-\tau \right)f\left(\tau \right)d\tau }$$ (1)where $\tau$ is the time-shift factor, and $\otimes$ represents the convolution operation. After obtaining $h\left(t\right)$ and $y\left(t\right)$ , $f\left(t\right)$ is calculated by the deconvolution of equation (1). Therefore, the force identification equation of the linear structure can be established by the deconvolution of equation (1), and the transform matrix is formed by the unit impulse response functions.

However, for a structure with the local nonlinearity studied in this article, the structural output response is not only generated by the excitation force but is also affected by the nonlinear force produced by nonlinear factors. ³⁵ Additionally, the actual frequency response function (FRF) is nonlinear FRF. As a result, for the local nonlinear substructure, $y\left(t\right)$ cannot be expressed as equation (1). This illustrates that the corresponding force identification equation cannot be obtained as a linear structure.

A local nonlinear substructure means that the number of nonlinear elements to be analyzed is far lower than the number of degrees of freedom of the entire structure. It also means that different from the structural nonlinearity caused by material or large deformation, local nonlinearity (which can be described by the lumped spring-damping element) only acts on the local position of the structure. Therefore, based on the output feedback principle, ³⁶ the nonlinear force generated by the local nonlinear factor and excitation force are regarded as the external force acting on the underlying linear part of the nonlinear structure. Then, for a structure with the local nonlinearity, the differential equation of motion can be rewritten as equations (2) and (3), which form the nonlinear state-space model^31,32 under defining state variable $\mathbf{z}\left(t\right)$ and extended force vector ${\mathbf{f}}_{ex}\left(t\right)$ $$\dot{\mathbf{z}}\left(t\right)={\mathbf{A}}_c\mathbf{z}\left(t\right)+{\mathbf{B}}_c{\mathbf{f}}_{ex}\left(t\right)={\mathbf{A}}_c\mathbf{z}\left(t\right)+{\mathbf{B}}_c\left[\begin{array}{c}f\left(t\right)\\ {\mathbf{g}}_{nl}\left(t\right)\end{array}\right]$$ (2) $$\mathbf{y}\left(t\right)={\mathbf{C}}_c\mathbf{z}\left(t\right)+{\mathbf{D}}_c{\mathbf{f}}_{ex}\left(t\right)$$ (3)in which $$\{\begin{array}{l}\mathbf{z}\left(t\right)={\left\{\begin{array}{cc}\mathbf{x}\left(t\right)& \dot{\mathbf{x}}\left(t\right)\end{array}\right\}}^{\text{T}}\\ {\mathbf{A}}_c=\left[\begin{array}{cc}{\mathbf{0}}_{n\times n}& {\mathbf{I}}_{n\times n}\\ -{\mathbf{M}}^{-1}\mathbf{K}& -{\mathbf{M}}^{-1}{\mathbf{C}}_{\alpha }\end{array}\right]\\ {\mathbf{B}}_c=\left[\begin{array}{cccc}\begin{array}{c}{\mathbf{0}}_{n\times 1}\\ {\mathbf{M}}^{-1}{\mathbf{L}}_f\end{array}& \begin{array}{c}{\mathbf{0}}_{n\times 1}\\ {\mathbf{M}}^{-1}{\mu }_1{\mathbf{L}}_{nl\_1}\end{array}& \begin{array}{c}\cdots \\ \cdots \end{array}& \begin{array}{c}{\mathbf{0}}_{n\times 1}\\ {\mathbf{M}}^{-1}{\mu }_p{\mathbf{L}}_{nl\_p}\end{array}\end{array}\right]\\ {\mathbf{g}}_{nl}\left(t\right)={\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}-g_1\left(t\right)& -g_2\left(t\right)\end{array}& \cdots \end{array}& -g_p\left(t\right)\end{array}\right]}^{\text{T}}\\ f_{nl\_i}\left(t\right)={\mu }_i{g}_i\left(t\right),\hspace{1em}i=\mathrm{1,2},\dots ,p\end{array}$$ (4)where ${\mathbf{A}}_c$ , ${\mathbf{B}}_c$ , ${\mathbf{C}}_c$ , and ${\mathbf{D}}_c$ denote the coefficient matrices of the continuous-time model. ${\text{C}}_{\alpha }$ , $\text{K}$ , and $\text{M}$ represent damping, stiffness, and mass matrices, respectively. $p$ is the total number of nonlinear factors. $f_{nl\_i}\left(t\right)$ represents the nonlinear force caused by the i-th nonlinear factor. ${\mu }_i$ and $g_i\left(t\right)$ are the corresponding nonlinear parameter and nonlinear describing function, respectively. ${\text{L}}_f$ and ${\text{L}}_{nl\_i}$ (composed of 0, -1 and 1) are the position vectors of excitation force $f\left(t\right)$ and nonlinear force $f_{nl\_i}\left(t\right)$ (or nonlinear describing function $g_i\left(t\right)$ ). For the nonlinear state-space model, extended force ${\mathbf{f}}_{ex}\left(t\right)$ including ${\mathbf{g}}_{nl}\left(t\right)$ and $f\left(t\right)$ is the input vector. Additionally, coefficient matrices ${\mathbf{C}}_c$ and ${\mathbf{D}}_c$ in equation (3) are the output influence matrices related to the type of output or measured response $\mathbf{y}\left(t\right)$ . For the stiffness nonlinearity studied in this article, the corresponding $f\left(t\right)$ or $g\left(t\right)$ is a function of the displacement response at nonlinear positions, so the displacement response is used to form the output response $\mathbf{y}\left(t\right)$ . Output influence matrices are set as ${\mathbf{C}}_c=\left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\end{array}\right]$ and ${\mathbf{D}}_c=\mathbf{0}$ , where $\mathbf{I}$ denotes the identity matrix.

Since vibration signals measured in engineering practice are discrete data corresponding to sampling time points, the discrete-time model shown in equations (5) and (6) needs to be developed from the corresponding continuous-time model shown in equations (2) and (3) ³⁷ $$\mathbf{z}\left(\left(k_t+1\right)\mathrm{\Delta }t\right)={\mathbf{A}}_d\mathbf{z}\left(k_t\mathrm{\Delta }t\right)+{\mathbf{B}}_d{\mathbf{f}}_{ex}\left(k_t\mathrm{\Delta }t\right)$$ (5) $$\mathbf{y}\left(k_t\mathrm{\Delta }t\right)={\mathbf{C}}_d\mathbf{z}\left(k_t\mathrm{\Delta }t\right)+{\mathbf{D}}_d{\mathbf{f}}_{ex}\left(k_t\mathrm{\Delta }t\right)$$ (6)where $k_t$ represents a non-negative integer, and $\mathrm{\Delta }t$ denotes the sampling interval. The corresponding conversion relationship is shown as $$\{\begin{array}{l}{\mathbf{A}}_d=e^{{\mathbf{A}}_c{k}_t\mathrm{\Delta }t};\hspace{1em}{\mathbf{C}}_d={\mathbf{C}}_c\\ {\mathbf{B}}_d=\left(e^{{\mathbf{A}}_c{k}_t\mathrm{\Delta }t}-\mathbf{I}\right){\mathbf{A}}_c^{-1}{\mathbf{B}}_c;\hspace{1em}{\mathbf{D}}_d={\mathbf{D}}_c\end{array}$$ (7)

For the mechanical system, an output response is generated only when there is an input signal. Therefore, zero initial condition $\mathbf{z}\left(0\right)=0$ is considered for such a system. Then, substituting the state variable $\mathbf{z}\left(k_t\mathrm{\Delta }t\right)$ ( $k_t=\mathrm{1,2},\dots ,N$ , $N$ is the total number of sampling points) described by equations (5) into (6), the output response vector is rewritten as $$\left[\begin{array}{c}\mathbf{y}\left(0\right)\\ \mathbf{y}\left(\mathrm{\Delta }t\right)\\ ⋮\\ \mathbf{y}\left(\left(N-1\right)\mathrm{\Delta }t\right)\end{array}\right]=\left[\begin{array}{cccc}{\mathbf{D}}_d& \mathbf{0}& \cdots & \mathbf{0}\\ {\mathbf{C}}_d{\mathbf{B}}_d& {\mathbf{D}}_d& \cdots & \mathbf{0}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{C}}_d{\mathbf{A}}_d^{N-2}{\mathbf{B}}_d& \cdots & {\mathbf{C}}_d{\mathbf{B}}_d& {\mathbf{D}}_d\end{array}\right]\left[\begin{array}{c}{\mathbf{f}}_{ex}\left(0\right)\\ {\mathbf{f}}_{ex}\left(\mathrm{\Delta }t\right)\\ ⋮\\ {\mathbf{f}}_{ex}\left(\left(N-1\right)\mathrm{\Delta }t\right)\end{array}\right]$$ (8)

The corresponding compact matrix-vector form of equation (8) is given by $$\mathbf{Y}=\mathbf{T}{\mathbf{F}}_{ex}$$ (9)

As shown in equation (9), output vector $\mathbf{Y}$ is directly related to extended force vector ${\mathbf{F}}_{ex}$ by transfer matrix $\mathbf{T}$ . Once $\mathbf{Y}$ and $\mathbf{T}$ have been obtained, unknown ${\mathbf{F}}_{ex}$ is obtained by solving equation (9). Then, unknown excitation force $f\left(t\right)$ is extracted from obtained vector ${\mathbf{F}}_{ex}$ . Therefore, in this study, equation (9) is viewed as the force identification equation of nonlinear structure in time domain.

The dynamic response can be measured to form the response vector $\mathbf{Y}$ . The transfer matrix $\mathbf{T}$ in equation (8) consists of these coefficient matrices, which are composed of the identity matrix, zero matrix, and structural characteristic matrices in equation (4). Therefore, the transfer matrix $\mathbf{T}$ is constant for different excitation cases on the same structure. In this study, these coefficient matrices can be identified by NSI algorithm^31,32 to form the transfer matrix $\mathbf{T}$ as equation (8).

Nonlinear subspace identification algorithm evolved from the classic linear subspace identification (LSI) algorithm ³⁸ and is a time-domain nonlinear system identification method. Based on equations (5) and (6), once output and input vectors have been obtained, NSI algorithm mainly containing row space transformation technology and matrix projection technology is used to calculate those coefficient matrices ( ${\mathbf{A}}_d$ , ${\mathbf{B}}_d$ , ${\mathbf{C}}_d$ , and ${\mathbf{D}}_d$ ). Details of the calculation procedure are given in Refs. 31 and 32. Then, identified coefficient matrices are applied to form the transfer matrix $\mathbf{T}$ . Furthermore, these discrete-time coefficient matrices can be applied to obtain the corresponding continuous-time coefficient matrices by equation (7). Then, the extended FRF ${\mathbf{H}}_{ex}\left(\omega \right)$ can be calculated by $$\begin{array}{l}{\mathbf{H}}_{\text{ex}}\left(\omega \right)={\mathbf{C}}_c{\left(\text{i}\omega \mathbf{I}-{\mathbf{A}}_c\right)}^{-1}{\mathbf{B}}_c+{\mathbf{D}}_c\\ =\mathbf{I}{\left(\mathbf{K}+\text{i}\omega {\mathbf{C}}_{\alpha }-{\omega }^2\mathbf{M}\right)}^{-1}\left[\begin{array}{cccc}{\mathbf{L}}_f& {\mu }_1{\mathbf{L}}_{nl\_1}& \cdots & {\mu }_p{\mathbf{L}}_{nl\_p}\end{array}\right]\hspace{0.17em}\\ ={\mathbf{H}}_L\left(\omega \right)\left[\begin{array}{cccc}{\mathbf{L}}_f& {\mu }_1{\mathbf{L}}_{nl\_1}& \cdots & {\mu }_p{\mathbf{L}}_{nl\_p}\end{array}\right]\hspace{0.17em}\end{array}$$ (10)where ${\mathbf{H}}_L\left(\omega \right)$ is the linear FRF corresponding to the underlying linear structure. This means that, in addition to forming a transfer matrix $\mathbf{T}$ , NSI algorithm can be introduced to calculate the linear FRF ${\mathbf{H}}_L\left(\omega \right)$ and nonlinear parameter $\mu$ by equation (10).

Regularization Method for Inverse Problem

Once transfer matrix $\mathbf{T}$ and output response vector $\mathbf{Y}$ are obtained, extended force vector ${\mathbf{F}}_{ex}$ can be calculated by solving equation (9) to eventually obtain the unknown excitation force $f\left(t\right)$ . According to the above information, it seems easy to conduct force identification of the nonlinear structure by the direct matrix inversion method. However, for the mechanical structure, the transfer matrix $\mathbf{T}$ is ill-conditioned because its condition number is always large. The direct inversion operation of ill-conditioned matrix $\mathbf{T}$ greatly amplifies the errors in the structural dynamic response and eventually obtains a worthless solution. Therefore, some classic regularization methods have been developed to solve this problem.

The measured dynamic response is inevitably polluted by interference information. Additionally, the estimated coefficient matrices are not absolutely accurate because they are identified on the basis of contaminated dynamic response, and NSI algorithm also produces the systematic errors. This indicates errors in both the transfer matrix $\mathbf{T}$ and response vector $\mathbf{Y}$ . Therefore, the TTLS^33,34 method is utilized to solve the force identification equation of nonlinear structures in this study.

Similar to TSVD method, TTLS method aims at neglecting small singular values. $\mathbf{T}$ and $\mathbf{Y}$ polluted by errors are combined to form an augmented matrix $\left(\mathbf{T},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathbf{Y}\right)$ , and the corresponding singular value decomposition is given by $$\left(\mathbf{T},\mathbf{Y}\right)=\mathbf{U}\mathbf{S}{\mathbf{V}}^{\text{T}}={\displaystyle {\sum }_{i=1}^{N+1}{\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^{\text{T}}}$$ (11)where $\mathbf{U}=\left[{\mathbf{u}}_1,{\mathbf{u}}_2,\dots ,{\mathbf{u}}_{N+1}\right]$ and $\mathbf{V}=\left[{\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_{N+1}\right]$ are the left singular matrix and right singular matrix of the augmented matrix, respectively. Element ${\sigma }_i$ in diagonal matrix $\mathbf{S}$ denotes the singular value.

Then, based on determining an appropriate truncation order (or threshold) $\eta$ , matrix $\mathbf{V}\in {ℝ}^{\left(N+1\right)\times \left(N+1\right)}$ can be partitioned as $$\mathbf{V}=\begin{array}{cc}\begin{array}{c}\\ N\\ 1\end{array}& \begin{array}{c}\begin{array}{cc}\eta & N-\eta +1\end{array}\\ \left(\begin{array}{cc}{\mathbf{V}}_{11}& {\mathbf{V}}_{12}\\ {\mathbf{V}}_{21}& {\mathbf{V}}_{22}\end{array}\right)\end{array}\end{array},\hspace{1em}\mathrm{subject}\hspace{0.17em}\mathrm{to}\{\begin{array}{c}{\mathbf{V}}_{11}\in {ℝ}^{N\times \eta }\\ {\mathbf{V}}_{22}\in {ℝ}^{1\times \left(N+1-\eta \right)}\end{array}$$ (12)

Finally, based on equation (12), TTLS solution of the force identification equation of the nonlinear structure is $${\mathbf{F}}_{ex\_TTLS}=-{\mathbf{V}}_{12}{\mathbf{V}}_{22}^{†}=-\frac{{\mathbf{V}}_{12}{\mathbf{V}}_{22}^{\text{T}}}{\Vert {\mathbf{V}}_{22}{\Vert }_2^2}$$ (13)where the superscript ${·}^{†}$ denotes the pseudo-inverse operation. As shown in equation (2), after obtaining the extended force vector ${\mathbf{F}}_{ex\_TTLS}$ , the excitation force acting on the nonlinear structure can eventually be identified.

It is noted from the above discussions that truncation order $\eta$ is regarded as the regularization parameter. This is a key parameter to improve the stability of solutions of ill-posed force identification problems. Many excellent approaches have been developed to obtain this parameter.^4,5,7 In this study, GCV criterion is applied to obtain an optimal truncation order $\eta$ of TTLS method. A description of the process can be found in Ref. 39.

In summary, the proposed excitation force identification strategy of the nonlinear structure contains two stages: an off-line preprocessing stage and a force identifying stage. The known excitation force generated by the vibration exciter is imposed on the nonlinear structure in the off-line preprocessing stage. Then, under the measured dynamic response and extended force (nonlinear describing function can be calculated by the dynamic response), the transfer matrix $\mathbf{T}$ can be formed by these coefficient matrices which are calculated by NSI algorithm. For the force identifying stage, a nonlinear structure operates under the actual (unknown) excitation force. TTLS method in conjunction with GCV criterion is applied to solve equation (9) and eventually identify the excitation force of the nonlinear structure. Furthermore, transfer matrix $\mathbf{T}$ should be a positive definite matrix or over definite matrix to facilitate solving equation (9). Therefore, the number of extended forces containing nonlinear describing functions and excitation force needs to be no more than the number of response measurement points. Additionally, NSI algorithm was originally used to identify nonlinear parameter ${\mu }_i$ on the basis of equation (10). Therefore, as shown in equations (2) and (4), ${\mu }_i$ is contained in the coefficient matrix ${\mathbf{B}}_c$ , and nonlinear describing function ${\mathbf{g}}_{nl}\left(t\right)$ is contained in extended force ${\mathbf{f}}_{ex}\left(t\right)$ . To identify the excitation force that is the focus of this study, ${\mu }_i$ can be assumed to have been identified by various nonlinear parameter identification algorithms. ⁴⁰ Then, ${\mu }_i$ can be moved into ${\mathbf{f}}_{ex}\left(t\right)$ and multiplied by the nonlinear description function to obtain the corresponding nonlinear force $f_{nl\_i}\left(t\right)={\mu }_i{g}_i\left(t\right)$ . This means that the nonlinear force vector ${\mathbf{f}}_{nl}\left(t\right)$ can replace ${\mathbf{g}}_{nl}\left(t\right)$ to form a new extended force vector ${\tilde{\mathbf{F}}}_{\mathrm{e}\mathrm{x}}$ in equation (9). Therefore, the same strategy can be used to calculate a new vector ${\tilde{\mathbf{F}}}_{\mathrm{e}\mathrm{x}}$ , and nonlinear force can also be extracted from ${\tilde{\mathbf{F}}}_{\mathrm{e}\mathrm{x}}$ similar to extracting the excitation force. This work is discussed in the following section.

Numerical and Experimental Studies

Two numerical simulation cases and one experimental case are utilized to verify the performance of the force identification strategy of the nonlinear structure in this section. Case 1 focuses on showing the feasibility of the proposed strategy and discussing the influence of nonlinear factors on the force identification result. Case 2 demonstrates the performance of the proposed strategy on multiple nonlinearities. For Case 3, the effectiveness of proposed strategy under actual experimental conditions is validated on a cantilever test bed with clearance nonlinearity.

In subsequent case studies, the force identification error ${\epsilon }_f$ defined by equation (14) is used to evaluate the closeness between the identification result ${\mathbf{f}}_{\text{iden}}$ and true value ${\mathbf{f}}_{\text{true}}$ $${\epsilon }_f=\frac{{\Vert {\mathbf{f}}_{\text{true}}-{\mathbf{f}}_{\text{iden}}\Vert }_2}{{\Vert {\mathbf{f}}_{\text{true}}\Vert }_2}\times 100\%$$ (14)

Furthermore, considering that the measured actual response is inevitably polluted by the interference information, the simulated noise signal ${\mathbf{e}}_n$ is added to the output response $\mathbf{y}$ for Case 1 and Case 2, and the contaminated output response ${\mathbf{y}}_n$ is expressed as $${\mathbf{y}}_n=\mathbf{y}+{\mathbf{e}}_n=\mathbf{y}+L_n\cdot \text{std}\left(\mathbf{y}\right)\cdot {\mathbf{\delta }}_n$$ (15)where $L_n$ denotes the noise level, and it can be mutually converted with the signal-to-noise ratio (SNR) commonly used in the signal processing field. ${\mathbf{\delta }}_n\in {ℝ}^{N\times 1}$ represents a pseudo-random number that obeys a uniform distribution in the interval (−1, 1).

Case 1: Five-DOF Structure with Cubic Stiffness Nonlinearity

As shown in Figure 1, a five-DOF simulation model with cubic stiffness nonlinearity is utilized to validate the feasibility of the force identification strategy of nonlinear structures. Excitation force acts at the third DOF. $m_1$ , $k_1$ , and $c_1$ in the figure are the mass, stiffness, and damping, respectively. $k_{cb}$ is a nonlinear parameter that denotes cubic stiffness. These structural parameters of the simulation model are set as $$\{\begin{array}{l}{m}_1=m_2=m_3=m_4=m_5=1\hspace{0.17em}\text{kg}\\ c_1=c_3=c_5=c_7=45\hspace{0.17em}\text{N}\cdot \text{s}/\text{m}\\ c_2=c_4=c_6=c_8=15\hspace{0.17em}\text{N}\cdot \text{s}/\text{m}\\ k_1=k_3=k_5=k_7=1.8\times {10}^4\left(\text{N}/\text{m}\right)\\ k_2=k_4=k_6=k_8=1.2\times {10}^4\left(\text{N}/\text{m}\right)\\ k_{cb}=2\times {10}^8{\left(\text{N}/\text{m}\right)}^3\end{array}$$ (16)

Figure 1.

Five-degree-of-freedom simulation model with cubic stiffness nonlinearity.

The cubic stiffness nonlinearity shown in Figure 1 is located between the fourth DOF and ground, and the corresponding nonlinear force is expressed as $$f_{nl\_cb}\left(t\right)=k_{cb}{g}_{cb}\left(t\right)=k_{cb}{x}_4^3\left(t\right)$$ (17)where $x_4\left(t\right)$ denotes the displacement response at the fourth DOF. The simulated output response is calculated by the fourth-order Runge–Kutta algorithm in MATLAB, which provides a data basis for force identification. The sampling frequency is 1024 Hz, and the acquisition duration is 2 s. The calculated output responses in the off-line preprocessing stage and force identifying stage are added by additive noise interference, and the noise level $L_n$ is 0.15.

The simulation case contains an excitation force and a nonlinear factor. Therefore, it is necessary to set two response measurement points (two reference locations). To conveniently calculate the nonlinear describing function of cubic stiffness nonlinearity by equation (17), the displacement responses at the second DOF and fourth DOF are used to conduct the force identification. Furthermore, according to the proposed force identification strategy of nonlinear structures, the known Gauss white noise signal is applied to conduct vibration testing in the off-line preprocessing stage. Gaussian white noise is produced by the “band-limited white noise” module in SIMULINK module of MATLAB. The “sample time” is set to 1/1024 s, and “noise power” is set to 0.1. Then, coefficient matrices of the nonlinear state-space model are calculated by NSI algorithm and are applied to obtain the transfer matrix. Based on equation (10), these coefficient matrices can also be used to estimate the extended FRF and eventually obtain the linear FRF and nonlinear parameter. As shown in Figure 2, the linear FRF between the excitation location (third DOF) and nonlinear location (fourth DOF) calculated by NSI algorithm is close to the corresponding theoretical value. Under vibration data of the underlying linear structure (removing nonlinear factors), the theoretical linear FRF is calculated by the classical H1 estimation method. ⁴¹ Additionally, as shown in Figure 3, the mean value of the identified cubic stiffness at 0–50 Hz is $2.037\times {10}^8{\left(\text{N}/\text{m}\right)}^3$ , which is close to the true value of $2\times {10}^8{\left(\text{N}/\text{m}\right)}^3$ . The above discussions indicate that these coefficient matrices identified by NSI algorithm are reliable and can be used to obtain the transfer matrix as equation (8).

Figure 2.

Linear frequency response function between excitation location and nonlinear location.

Figure 3.

Identification result of cubic stiffness by nonlinear subspace identification algorithm.

In the force identification stage, the excitation force is set as $f\left(t\right)=20sin\left(2\pi \times 15t\right)\times \left[1-\cos \left(2\pi t\right)\right]$ . As long as the excitation location and response measurement location remain unchanged, the transfer matrix is constant for different excitation cases on the same structure. The condition number of the transfer matrix is $1.43\times {10}^{16}$ . This shows that the transfer matrix is an ill-conditioned matrix. Therefore, for the above nonlinear simulation model, the force identification problem is ill-posed, and it is solved by TTLS method combined with GCV criterion in this study. As shown in Figure 4, the optimal truncation order $\eta$ is 2087. Then, Figure 5 shows the identification result of the excitation force. The red dotted line represents the true excitation force calculated by the corresponding function expression. The identified excitation force curve is close to the true value curve, and the excitation force identification error is 5.56%. This illustrates that the proposed identification strategy is suitable for identifying the excitation force of nonlinear structures.

Figure 4.

Generalized cross-validation criterion for choosing truncation order.

Figure 5.

Identification result of excitation force.

When the above nonlinear simulation model is linearized by ignoring the existing cubic stiffness nonlinearity, the identified excitation force is shown in Figure 6. In the identification process, the force identification equation of a linearized structure is established on the basis of the corresponding linear state-space model. Nonlinear parameters are not included in matrix ${\mathbf{B}}_c$ , and the extended force vector contains only the excitation force.^11,12 Instead of NSI method, the LSI method ³⁸ is applied to estimate new coefficient matrices to obtain a new transfer matrix. Then, the solving method is also TTLS combined with GCV criterion. Compared to Figure 5, Figure 6 shows an obvious difference in the region with the largest amplitude. The identification error is 11.97%. This illustrates that to obtain an accurate excitation force of nonlinear structures, existing nonlinear factors cannot be ignored.

Figure 6.

Identification result of excitation force when ignoring nonlinear factor.

It is worth noting that in addition to the excitation force, a nonlinear describing function can also be abstracted from the calculated extended force vector. However, as shown in Figure 7, the identification result of the nonlinear describing function seriously deviates from its true value. This is because the magnitude order of nonlinear describing function (approximately ${10}^{-7}$ ) is much smaller than that of the excitation force (approximately $10$ ). The calculation errors corresponding to the above two parts are coupled. This can submerge the identification result of the nonlinear describing function. As discussed in section 3, nonlinear parameter ${\mu }_i$ contained in matrix ${\mathbf{B}}_c$ can be moved into ${\mathbf{f}}_{ex}\left(t\right)$ . This means that instead of the nonlinear describing function, nonlinear force can be contained in the extended force vector. The magnitude order of the nonlinear force is close to that of the excitation force. Figure 8 shows the identification result of nonlinear force with the same simulation case. When ill-posed equation (9) is solved by TTLS method, the solutions are essentially obtained under the energy minimization principle. Therefore, the solutions are always nonzero over the entire time range,^12,13 and there is a difference between the identification result and true value in the area where the true value is equal to 0, as shown in Figure 8. However, in other areas, the identification result is close to the true value. This further illustrates that the proposed identification strategy is feasible.

Figure 7.

Identification result of nonlinear describing function.

Figure 8.

Identification result under nonlinear force considered in extended force vector.

Case 2: Cantilever Model with Two Nonlinearities

In this section, a cantilever model with two square nonlinearities shown in Figure 9 is built to verify the performance of the identification strategy on a continuous structure with multiple nonlinearities. In the figure, two kinds of square nonlinearities are located on the simulated cantilever beam. The corresponding nonlinear forces are given by $$f_{nl\_sn1}\left(t\right)=k_{sn1}{g}_{sn1}\left(t\right)=k_{sn1}{x}^2\left(l_{sn1},t\right)$$ (18) $$f_{nl\_sn2}\left(t\right)=k_{sn2}{g}_{sn2}\left(t\right)=k_{sn1}x\left(l_{sn2},t\right)\cdot \left|x\left(l_{sn2},t\right)\right|$$ (19)where $x\left(l_{sn1},t\right)$ and $x\left(l_{sn2},t\right)$ are the displacement responses at two nonlinear positions. $k_{sn1}$ and $k_{sn2}$ are the nonlinear stiffness parameters of the above two square nonlinearities.

Figure 9.

Simulated cantilever beam with two square nonlinearities.

In the simulation calculation, the Young’s modulus, density, cross-sectional area, cross-sectional moment of inertia, and length of the cantilever model are $2.1\times {10}^{11}\text{Pa}$ , $7.85\times {10}^3{\left(\text{kg}/\text{m}\right)}^3$ , $2\times {10}^{-4}{\text{m}}^2$ , $1.67\times {10}^{-9}{\text{m}}^4$ , and 1 m, respectively. The nonlinear stiffness parameters are set as $k_{sn1}=2\times {10}^7{\left(\text{N}/\text{m}\right)}^2$ and $k_{sn2}=4\times {10}^8{\text{N}/\text{m}}^2$ . The acting positions of excitation force and two nonlinear factors are set as $l_f=0.2\hspace{0.17em}\text{m}$ , $l_{sn1}=0.6\hspace{0.17em}\text{m}$ , and $l_{sn2}=0.85\hspace{0.17em}\text{m}$ , respectively. The acquisition duration and sampling frequency are 1 s and 2048 Hz, respectively. Additionally, the noise level $L_n$ is set as 0.20.

To conduct force identification of the above simulation model, three response measurement points are set at the excitation location and two nonlinear locations. Additionally, in the off-line preprocessing stage, Gaussian white noise is produced by “band-limited white noise” module in the SIMULINK module of MATLAB. “Sample time” is set to 1/2048 s, and “noise power” is set to 0.1. After forming the transfer matrix in the off-line preprocessing stage, excitation force $f\left(t\right)=10sin\left(2\pi \times 20t\right)\times \left[2-\cos \left(2\pi \times 5t\right)\right]$ is used to conduct vibration testing on the simulation model in the force identifying stage. Figure 10(a) shows the excitation force identification result by TTLS method, and the identification error is 9.41%. The TSVD method is also a regularization method for solving ill-posed identification equations by neglecting the small singular values of the transfer matrix. When equation (9) is solved by TSVD method, the identified excitation force is shown in Figure 10(b). There is an obvious difference between the true force curve and identified excitation force curve, and the identification error is 16.46%. This occurs because compared with TTLS method, TSVD method only considers errors in the response vector $\mathbf{Y}$ . However, in the off-line preprocessing stage, the dynamic response for obtaining the transfer matrix is also inevitably polluted by noise interference. Furthermore, there is currently no algorithm that can identify nonlinear systems with absolute accuracy. This means that errors should be considered in both transfer matrix $\mathbf{T}$ and response vector $\mathbf{Y}$ . Consequently, Figure 10 shows that the identification error by TTLS method is smaller than that by TSVD method. Additionally, when two nonlinear forces are considered in the extended force vector, nonlinear force identification results by TTLS method are shown in Figure 11. Similar to the excitation force identification results, the identified nonlinear force curves are close to the true value curves. The identification errors of the two nonlinear forces are 9.89% and 12.65%, respectively. The above discussions demonstrate that the proposed identification strategy is suitable for a continuous structure with nonlinearities.

Figure 10.

Excitation force identification results: (a) truncated total least squares method; (b) truncated singular value decomposition method.

Figure 11.

Nonlinear force identification results: (a) $f_{nl\_sn1}\left(t\right)$ ; (b) $f_{nl\_sn2}\left(t\right)$ .

Case 3: Cantilever Test Bed with Clearance Nonlinearity

As shown in Figure 12, the experimental verification of the proposed force identification strategy of nonlinear structures is conducted on a cantilever test bed with clearance nonlinearity. The test bed includes a clearance device and cantilever beam. The clearance device is composed of two small cantilever beams and a fixed base. Additionally, as shown in the enlarged part of the figure, clearance nonlinearity in this test bed means that a gap exists between the contact head and cantilever. In the experiment, a vibration exciter is utilized to impose the excitation force acting on the test bed, and an impedance head which is a force sensor installed on the vibration exciter is used to measure the true excitation force. Eddy current sensors are installed at some measurement positions, and the measured displacement responses are used to identify the excitation force.

Figure 12.

Cantilever test bed with clearance nonlinearity.

The clearance nonlinear force caused by clearance nonlinearity is given by $$f_{nl\_c}\left(t\right)=k_c{g}_c\left(t\right)=k_c·\{\begin{array}{c}x(l_c,t)-d_{c}x(l_c,t)\ge d_c\\ 0-d_c<x(l_c,t)<d_c\\ x(l_c,t)+d_{c}x(l_c,t)\le -d_c\end{array}$$ (20)where $k_c$ and $d_c$ denote the clearance stiffness and clearance value, respectively. $x\left(l_c,t\right)$ represents the displacement response at the clearance location. For the proposed identification strategy, the nonlinear describing function should be obtained to form the extended force vector. Equation (20) shows that the clearance describing function $g_c\left(t\right)$ includes the displacement response $x\left(l_c,t\right)$ and clearance value $d_c$ . Therefore, $x\left(l_c,t\right)$ and $d_c$ should be obtained first. The dynamic measurement method and identification method of the clearance value were fully discussed in our previous studies. ³⁵

In the experiment, the excitation position and clearance nonlinearity are located at $l_f=0.1\hspace{0.17em}\text{m}$ and $l_c=0.575\hspace{0.17em}\text{m}$ . Two response measurement points are set at $l_{m1}=0.3\hspace{0.17em}\text{m}$ and $l_{m2}=0.575\hspace{0.17em}\text{m}$ . In the off-line preprocessing stage, Gaussian white noise is generated by the signal generator, the frequency band is 0–400 Hz, and the corresponding signal intensity or power is set to 512 W, wherein the amplitude of PSD (power spectral density) is set to 0.25. The signal is amplified by a power amplifier and then used to conduct vibration testing to obtain the transfer matrix. Additionally, as shown in Figure 13, the linear FRF ( ${\mathbf{H}}_{L21}$ and ${\mathbf{H}}_{L31}$ ) of the underlying linear part of nonlinear structures between the excitation position (point 1) and two measurement positions (point 2 and point 3) can also be obtained in this stage. Then, in the force identifying stage, the excitation force is set as a sinusoidal signal at 15 Hz. Given the displacement responses at two measurement points, Figure 14 shows the identification results of the excitation force by TTLS and TSVD methods, and the corresponding identification errors are 7.70% and 13.48%. This further illustrates that, in practice, it is difficult to obtain an absolutely accurate transfer matrix of nonlinear structures, and TLLS method is suitable for the force identification problem in which both the transfer matrix $\mathbf{T}$ and response vector $\mathbf{Y}$ contain errors.

Figure 13.

Linear frequency response function of underlying linear structure. (a) ${\mathbf{H}}_{L21}$ ; (b) ${\mathbf{H}}_{L31}$ .

Figure 14.

Excitation force identification results. (a) Truncated total least squares method; (b) truncated singular value decomposition method.

To further validate the performance of the proposed identification strategy, a new experiment is conducted by adjusting the excitation position and measurement position. The excitation position and clearance nonlinearity are located at $l_f=0.3\hspace{0.17em}\text{m}$ and $l_c=0.575\hspace{0.17em}\text{m}$ . Two response measurement points are set at $l_{m1}=0.1\hspace{0.17em}\text{m}$ and $l_{m2}=0.575\hspace{0.17em}\text{m}$ . The excitation force is set as the sinusoidal signal at 50 Hz. Figure 15 shows the identification result of the excitation force. It can be observed from the figure that the identified excitation force curve is close to the measured curve, and the identification error is 9.83%. In conclusion, the above identification results illustrate that the proposed identification strategy is also effective for identifying the excitation force of nonlinear structures in actual experiments. Furthermore, the test bed is simplified from the artillery barrel-cradle structure with clearance. ³⁵ For artillery, abnormalities in the firing force reduce the firing precision. The experimental results also demonstrate that the proposed identification strategy makes it possible to obtain the firing force acting on the artillery barrel to monitor the artillery firing performance.

Figure 15.

Excitation force identification results.

Conclusions

This study focused on identifying the excitation force of nonlinear structures, and a novel force identification strategy was proposed. The strategy contains two identification stages: an off-line preprocessing stage and a force identifying stage. For the identification strategy, based on a nonlinear state-space model, the force identification equation of the nonlinear structure was built. Its transfer matrix consists of coefficient matrices of the nonlinear state-space model. The input vector is an extended force vector composed of the nonlinear describing function and excitation force. The purpose of the off-line preprocessing stage is to apply NSI algorithm to obtain the transfer matrix. Then, in the force identifying stage, the TTLSs method in conjunction with GCV criterion is introduced to solve the ill-posed identification equation to eventually obtain the excitation force.

Two numerical simulation cases (five-DOF structure with cubic stiffness nonlinearity and cantilever model with two square nonlinearities) and one experiment case (cantilever test bed with clearance nonlinearity) were utilized to verify the effectiveness and feasibility of the proposed force identification strategy. The simulation and experimental results demonstrated that the proposed force identification strategy is suitable for identifying the external force of nonlinear structures and is not limited by the number and type of nonlinearities. It is worth noting that the main contribution here established a force identification equation of the nonlinear structure in time domain and solved the problem of obtaining the transfer matrix. Furthermore, the proposed identification strategy can extend the regularization methods that were previously used for a linear structure to a structure with the local nonlinearity. However, this paper does not propose a novel regularization method. Although identification results by TTLS method are acceptable, there are several problems. In particular, the singular value decomposition operation of the augmented matrix $\left(\mathbf{T},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathbf{Y}\right)$ greatly increases the computational cost. Some scholars have studied how to improve the performance of TTLS method.^42,43 In future work, we will also focus on the regularization method for ill-posed inverse problems of nonlinear structures.