Simultaneous identification of the load and unknown parameters of the structure based on the perturbation method

Abstract

For when a structure is both subjected to unknown loads and characterized by unknown parameters and structural parameter identification and load identification algorithms cannot be applied under these circumstances, this article proposes a new algorithm based on the perturbation method for the simultaneous identification of the load and unknown structural parameters using a few response points. The impulse response matrix is expanded with respect to the unknown parameters, and then the load term is combined with the unknown parameter perturbation term to generate the parameter–load term; finally, the load and unknown parameters can be simultaneously identified by iteration. The proposed algorithm is verified by an example of a 10-storey shear frame, and the effect of noise level, the number of response points, and the truncation ratio, which is a parameter introduced to improve the accuracy of the proposed algorithm, are studied. Moreover, the effect of the distribution of response points is discussed for another example of a simply supported beam, and the results show that when the response points are distributed over the full beam, the error is obviously smaller than when the response points are distributed over only part of the beam.

Introduction

Civil engineering structures have been subjected to environmental and human loads during many decades of operation, which may cause damage to the structures, thereby changing the structural parameters. The structural vibration response is a direct representation of the structural parameters. Many damage identification algorithms ¹ determine the position and size of the damage by estimating the changes in parameters. Therefore, it is a significant research direction to monitor the response of the structure and then obtain the structural parameters using the response to evaluate the current structural security situation.

Traditional structural parameter identification algorithms require the input and output of the structure to obtain the frequency response function in frequency domain or impulse response function in time domain to identify structural parameters, which is obviously difficult to realize in most cases. The structural parameter identification algorithms based on environmental excitation ² do not require structural input, for example, time series analysis, ³ natural excitation technique (NExT), ⁴ and stochastic subspace identification (SSI), ⁵ but the accuracy of structural parameter identification algorithms decreases, and even the algorithms may fail when the external load does not satisfy the Gaussian distribution assumption. For structural dynamic load identification algorithms, ⁶ it is necessary to know all the parameter information of the structure accurately in advance. However, some structural parameters cannot be easily obtained, especially after the structure is damaged.

Considering that a structure has both unknown loads and unknown parameters, if the unknown loads can be identified simultaneously with the unknown structural parameters, not only can the structural safety condition be evaluated by comparison of the identified structural parameters and the initial design parameters, but also the load source can be roughly determined according to the time history of the identified load. Meanwhile, more appropriate measures can be taken to control the load source or change the transmission path of the load to reduce the structural response, which can effectively lower the risk to the structure to ensure its long-term stable and safe operation.

In recent years, only the use of structural dynamic response to identify loads and unknown parameters has gradually developed. The most common algorithms are based on least-squares method^7–10 or Kalman filtering.^11–14 Moreover, Zhu and Law ¹⁵ identified the damage of a concrete bridge under moving load, and the interaction forces between the vehicle and the bridge are also identified simultaneously. A damage function is adopted to simulate the crack damage. Lu and Law ¹⁶ proposed a composite algorithm based on the sensitivities of response with respect to the load and unknown parameters. Zhang et al.^17,18 solved the problem of the composite identification of the load and parameters based on the virtual distortion method (VDM). Huang et al. ¹⁹ developed a new algorithm, referred to as the adaptive quadratic sum–squares error with unknown excitations (AQSSE-UI) by minimizing the quadratic sum–squares error between the measured response and the theoretical values. Law and Ding ²⁰ proposed two sub-structural identification methods with the structure divided into substructures and with one substructure assessed at one time. Sun and Betti ²¹ presented a hybrid heuristic optimization strategy based on the artificial bee colony (ABC) algorithm with a local search operator. Feng et al. ²² proposed a new method to identify structural parameters and vehicle dynamic axle loads of a vehicle–ridge interaction system in the framework of an iterative parametric optimization process. A Bayesian inference regularization is presented to solve the ill-posed least-squares problem. Pioldi and Rizzi ²³ proposed a new element-level system identification and input estimation technique named full dynamic compound inverse method (FDCIM). A statistical average technique, a modification process, and a parameter projection strategy are adopted at each stage to achieve stronger convergence. Jayalakshmi and Rao ²⁴ used a newly developed dynamic hybrid adaptive firefly algorithm (DHAFA) to solve the inverse problem associated with the system identification.

Based on the perturbation theory, this article proposes a new algorithm for simultaneous identification of the load and unknown parameters of the structure using a few response points. First, initial values are assumed for unknown parameters, and the convolution equation for solving the structural dynamic response is discretized. Next, the impulse response matrix is expanded with respect to unknown parameters at the assumed initial values. The partial derivative matrices of the impulse response matrix with respect to the unknown parameters are obtained by the difference method. The load term is combined with each unknown parameter perturbation term to generate the parameter–load term, and then the load term and the parameter–load term can be determined by inversion with Tikhonov regularization method ²⁵ and L-curve criterion. ²⁶ Each unknown parameter perturbation term is calculated using the ratio of the norm of the corresponding parameter–load term to that of the load term, and its positive and negative proprieties are determined by the method proposed in this article. Subsequently, all unknown parameter perturbation terms are added to initial values, and the resulting values are treated as the new initial values. Finally, the next iteration step proceeds until the identified results converge. In an example of a 10-storey shear frame structure, the effect of the number of response points, noise level, and truncation ratio on the identified results are studied. The proposed algorithm is also applied to the identification of another example of a simply supported beam under the action of a concentrated load, and the effect of the distribution of response points is also discussed.

Simultaneous identification of the load and unknown parameters

For a single-degree-of-freedom system with damping, the equation of motion can be written as

$$m\stackrel{··}{y}(t)+c\stackrel{·}{y}(t)+\mathrm{ky}(t)=f(t)$$

(1)

in which, m, c, and k are the mass, damping, and stiffness, respectively; $\stackrel{··}{y}(t)$ , $\stackrel{·}{y}(t)$ , and $y(t)$ are the acceleration, velocity, and displacement response, respectively; and $f(t)$ is the external load.

Suppose that both the initial velocity and the initial displacement are zero, then the convolution equation for solving the structural displacement response can be obtained by Duhamel integral

$$y(t)=\underset{0}{\overset{t}{\int }}{e}^{-\xi \omega (t-\tau )}\sin ({\omega }_d(t-\tau ))f(\tau )d\tau$$

(2)

in which $\xi =c/(2\sqrt{\mathrm{km}})$ is the damping ratio, $\omega =\sqrt{k/m}$ is the natural frequency, and ${\omega }_d=\omega \sqrt{1-{\xi }^2}$ is the damped natural frequency.

Let $h(t)=e^{-\xi \omega t}\sin ({\omega }_{d}t)$ , which is the impulse response function, and equation (2) can be discretized according to the time sampling points

$$\begin{array}{c}{\left\{\begin{array}{c}y(t_1)\hfill \\ \multicolumn{1}{c}{y(t_2)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{y(t_Q)}\end{array}\right\}}_{Q\times 1}={\left[\begin{array}{cccc}h(t_1)\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{h(t_2)}& h(t_1)\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{h(t_Q)}& h(t_{Q-1})\hfill & \cdots \hfill & h(t_1)\hfill \end{array}\right]}_{Q\times Q}\hfill \\ {\left\{\begin{array}{c}f(t_0)\hfill \\ \multicolumn{1}{c}{f(t_1)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{f(t_{Q-1})}\end{array}\right\}}_{Q\times 1}\mathrm{\Delta }t\hfill \end{array}$$

(3)

in which $y(t_k)$ , $f(t_k)$ , and $h(t_k)$ are the values of the displacement, load, and impulse response function at the sampling time point $t_k$ , respectively; $t_k=k\mathrm{\Delta }t$ is the kth sampling time point; and $\mathrm{\Delta }t$ is the sampling time interval.

Similarly, for a N-degrees-of-freedom system with proportional damping, the equation of motion can be written as

$$\mathit{M}\stackrel{··}{\mathit{y}}(t)+\mathit{C}\stackrel{·}{\mathit{y}}(t)+\mathit{Ky}(t)=\mathit{B}\mathit{f}(t)$$

(4)

in which ${\mathit{M}}_{N\times N}$ , ${\mathit{C}}_{N\times N}$ , and ${\mathit{K}}_{N\times N}$ are the mass, damping, and stiffness matrices, respectively; $\mathit{C}=\alpha \mathit{M}+\beta \mathit{K}$ , where $\alpha$ and $\beta$ are the mass matrix and stiffness matrix correlation coefficients of the damping matrix, respectively; $\stackrel{··}{\mathit{y}}(t{)}_{N\times 1}$ , $\stackrel{·}{\mathit{y}}(t{)}_{N\times 1}$ , and $\mathit{y}(t{)}_{N\times 1}$ are the acceleration, velocity, and displacement vectors, respectively; $\mathit{f}(t{)}_{N_f\times 1}$ is the load vector; $N_f$ is the number of external loads; and ${\mathit{B}}_{N\times N_f}$ is the distribution matrix associated with $\mathit{f}(t)$ .

Suppose that the initial velocity and the initial displacement vectors are zero. Using the orthogonality of the mass matrix, stiffness matrix, and damping matrix with respect to the mode shape matrix, the displacement of the selected response points calculated by the superposition method is discretized according to time sampling points as in equation (3)

$$\begin{array}{c}{\left\{\begin{array}{c}{\mathit{y}}_1\hfill \\ \multicolumn{1}{c}{{\mathit{y}}_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathit{y}}_{N_d}}\end{array}\right\}}_{(N_d\times Q)\times 1}={\left[\begin{array}{cccc}{\mathit{H}}_{11}\hfill & {\mathit{H}}_{12}\hfill & \cdots \hfill & {\mathit{H}}_{1N_f}\hfill \\ \multicolumn{1}{c}{{\mathit{H}}_{21}}& {\mathit{H}}_{22}\hfill & \cdots \hfill & {\mathit{H}}_{2N_f}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\mathit{H}}_{N_{d}1}}& {\mathit{H}}_{N_{d}2}\hfill & \cdots \hfill & {\mathit{H}}_{N_d{N}_f}\hfill \end{array}\right]}_{(N_d\times Q)\times (N_f\times Q)}\hfill \\ {\left\{\begin{array}{c}{\mathit{f}}_1\hfill \\ \multicolumn{1}{c}{{\mathit{f}}_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathit{f}}_{N_f}}\end{array}\right\}}_{(N_f\times Q)\times 1}\hfill \end{array}$$

(5)

in which ${\mathit{y}}_i$ $(i=1,2,\dots ,N_d)$ is the displacement of the ith response point, $N_d$ is the number of response points, ${\mathit{f}}_j$ $(j=1,2,\dots ,N_f)$ is the jth load, and ${\mathit{H}}_{\mathrm{ij}}$ $(i=1,2,\dots ,N_d,j=1,2,\dots ,N_f)$ is the discrete impulse response matrix with respect to the ith response point and the jth load.

It is assumed that the structure contains l unknown parameters, which can be expressed as

$$\mathit{b}=[\begin{array}{cccc}{b}_1\hfill & b_2\hfill & \cdots \hfill & b_l\hfill \end{array}{]}^T$$

(6)

in which the superscript “T” represents the matrix transposition.

Equation (5) can be simplified into the following form containing the unknown parameters

$$\mathit{y}=\mathit{H}(\mathit{b})\mathit{f}$$

(7)

According to the perturbation theory, first, the initial values are assumed for all unknown parameters, and then the impulse response matrix $\mathit{H}(\mathit{b})$ is used to perform the first-order Taylor expansion with respect to all unknown parameters at the assumed initial values

$$\mathit{H}(\mathit{b})={\mathit{H}}_0+\delta \mathit{H}={\mathit{H}}_0+\sum _{i=1}^l{\mathit{H}}_i·\delta b_i$$

(8)

in which ${\mathit{H}}_0$ is the impulse response matrix at the initial values of unknown parameters, ${\mathit{H}}_i$ $(i=1,2,\dots ,l)$ is the first-order partial derivative matrix of ${\mathit{H}}_0$ with respect to the ith unknown parameter $b_i$ , and $\delta b_i$ is the ith unknown parameter perturbation term.

Substituting equation (8) into equation (7), one can obtain

$$\mathit{y}=({\mathit{H}}_0+\sum _{i=1}^l{\mathit{H}}_i·\delta b_i)\mathit{f}={\mathit{H}}_0\mathit{f}+\sum _{i=1}^l{\mathit{H}}_i·(\delta b_i·\mathit{f})$$

(9)

Combining the ith unknown parameter perturbation term $\delta b_i$ and the load term $\mathit{f}$ , named the ith parameter–load term, that is, ${\mathit{f}}_i$ , equation (9) can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{\mathit{y}={\mathit{H}}_0\mathit{f}+\sum _{i=1}^l{\mathit{H}}_i·{\mathit{f}}_i}\\ =[\begin{array}{cccc}{\mathit{H}}_0\hfill & {\mathit{H}}_1\hfill & \cdots \hfill & {\mathit{H}}_l\hfill \end{array}]{\left[\begin{array}{cccc}\mathit{f}\hfill & {\mathit{f}}_1\hfill & \cdots \hfill & {\mathit{f}}_l\hfill \end{array}\right]}^T\hfill \end{array}$$

(10)

Solving the generalized inverse of equation (10), the load term and parameter–load term can be obtained

$${\left[\begin{array}{cccc}\mathit{f}\hfill & {\mathit{f}}_1\hfill & \cdots \hfill & {\mathit{f}}_l\hfill \end{array}\right]}^T=[\begin{array}{cccc}{\mathit{H}}_0\hfill & {\mathit{H}}_1\hfill & \cdots \hfill & {\mathit{H}}_l\hfill \end{array}{]}^{+}\mathit{y}$$

(11)

in which the superscript “+” represents the generalized inverse.

Then, the ith unknown parameter perturbation term can be determined using the ratio of the norm of ${\mathit{f}}_i$ to that of $\mathit{f}$

$$\left|\delta b_i\right|=\frac{\mathrm{norm}({\mathit{f}}_i)}{\mathrm{norm}(\mathit{f})}(i=1,2,\dots ,l)$$

(12)

in which |·| represents the absolute value, and $\mathrm{norm}(·)$ represents the vector norm, which can be the 1-norm or 2-norm.

The absolute value of an unknown parameter perturbation term is given for the reason that the ratio of the norm of ${\mathit{f}}_i$ $(i=1,2,\dots ,l)$ to that of $\mathit{f}$ is always positive, whereas $\delta b_i$ is not, because the actual value of an unknown parameter may be greater than or less than its initial value; thus, the positive and negative properties of the unknown parameter perturbation terms obtained by equation (12) must be determined.

Suppose that the load term $\mathit{f}={\left\{\begin{array}{cccc}{f}_1\hfill & f_2\hfill & \cdots \hfill & f_Q\hfill \end{array}\right\}}^T$ and the ith parameter–load term ${\mathit{f}}_i={\left\{\begin{array}{cccc}{f}_{i1}\hfill & f_{i2}\hfill & \cdots \hfill & f_{\mathrm{iQ}}\hfill \end{array}\right\}}^T$ are obtained one time; first, the corresponding elements of ${\mathit{f}}_i$ and $\mathit{f}$ are divided one by one to obtain a new term named $\mathit{r}{\mathit{f}}_i$

$$\mathit{r}{\mathit{f}}_i={[\frac{f_{i1}}{f_1}\frac{f_{i2}}{f_2}\cdots \frac{f_{\mathrm{iQ}}}{f_Q}]}^T(i=1,2,\dots ,l)$$

(13)

Next, the elements in $\mathit{r}{\mathit{f}}_i$ are accumulated, and the positive and negative properties of the unknown parameter perturbation term $\delta b_i$ can be determined by the positive and negative properties of the cumulative result of $\mathit{r}{\mathit{f}}_i$ . However, because it is inevitable that the values of some elements in $\mathit{f}$ are very close to zero, the absolute values of some elements in $\mathit{r}{\mathit{f}}_i$ will be too large, and these elements may affect the final cumulative result. Thus, the following method is used to solve this problem.

First, $\mathit{r}{\mathit{f}}_i$ is sorted by numerical value, and an equal number of larger and smaller elements are truncated. The ratio of the number of larger or smaller elements truncated to the length of $\mathit{r}{\mathit{f}}_i$ is defined as the truncation ratio. Next, the remaining elements are accumulated to determine the positive and negative properties of the unknown parameter perturbation term. A discussion on the truncation ratio is provided in the example. It should also be noted that it is not recommended to use the number of positive and negative elements in $\mathit{r}{\mathit{f}}_i$ to determine the positive and negative properties of the unknown parameter perturbation term because this method does not yield a good result.

After the perturbation terms are determined, they are added to the initial values of the unknown parameters, and then the results are taken as the initial values in the next iteration step, until the identified results converge.

In addition, for a multiple-degrees-of-freedom structure, it is difficult, even impossible to obtain the explicit expression of the impulse response function with respect to a certain structural parameter. Therefore, the difference method is used. Specifically, suppose that an unknown parameter has a small perturbation at its initial value; then, the discrete impulse response functions are calculated before and after the occurrence of the perturbation, and the difference between the two impulse response functions is treated as the partial derivative of the impulse response function with respect to this unknown parameter.

Introduction of the orthogonal polynomials

It can be seen from equation (7) that when the load and unknown parameters are simultaneously identified, the number of unknown terms becomes $(l+1)$ times that when only the load is identified. To ensure that the matrix $[\begin{array}{cccc}{\mathit{H}}_0\hfill & {\mathit{H}}_1\hfill & \cdots \hfill & {\mathit{H}}_l\hfill \end{array}]$ satisfies the reversible condition, more response points must be selected, a task that is generally difficult to achieve in engineering. Thus, Legendre polynomials are introduced in this section ²⁷

$$\begin{array}{c}{L}_1(t)=1\hfill \\ L_2(t)=\frac{2t}{T_0}-1\hfill \\ L_{v+1}(t)=\frac{2v-1}{v}·(\frac{2t}{T_0}-1)·L_v(t)\hfill \\ -\frac{v-1}{v}·L_{v-1}(t)(v=2,3,\dots ,V-1)\hfill \end{array}$$

(14)

in which $T_0$ represents that the definition interval of the Legendre polynomials is $[0,T_0]$ .

The load term $\mathit{f}$ and the parameter–load term ${\mathit{f}}_i$ are orthogonally expanded with respect to the Legendre polynomials and then discretized, similar to equation (3)

$$\begin{array}{c}\mathit{f}=\mathit{L}{\mathit{a}}_0\hfill \\ {\mathit{f}}_i=\mathit{L}{\mathit{a}}_i(i=1,2,\dots ,l)\hfill \end{array}$$

(15)

in which $\mathit{L}$ is the Legendre polynomial matrix; ${\mathit{a}}_0$ and ${\mathit{a}}_i$ $(i=1,2,\dots ,l)$ are the orthogonal expansion coefficients of the load term and the ith parameter–load term with respect to the Legendre polynomials, respectively.

Substituting equation (15) into equation (10), one can obtain

$$\begin{array}{c}\multicolumn{1}{c}{\mathit{y}=\left[\begin{array}{cccc}{\mathit{H}}_0\hfill & {\mathit{H}}_1\hfill & \cdots \hfill & {\mathit{H}}_l\hfill \end{array}\right]{\left[\begin{array}{cccc}\mathit{L}{\mathit{a}}_0\hfill & \mathit{L}{\mathit{a}}_1\hfill & \cdots \hfill & \mathit{L}{\mathit{a}}_l\hfill \end{array}\right]}^T}\\ =\left[\begin{array}{cccc}{\mathit{H}}_0\mathit{L}\hfill & {\mathit{H}}_1\mathit{L}\hfill & \cdots \hfill & {\mathit{H}}_l\mathit{L}\hfill \end{array}\right]{\left[\begin{array}{cccc}{\mathit{a}}_0\hfill & {\mathit{a}}_1\hfill & \cdots \hfill & {\mathit{a}}_l\hfill \end{array}\right]}^T\hfill \\ =\mathit{Ha}\hfill \end{array}$$

(16)

In general, it is usually sufficient to use a few orthogonal polynomials to accurately reconstruct the load, and thus, the number of unknown terms that must be identified is greatly reduced, enabling the identification of the load and unknown parameters using fewer response points. Therefore, it is difficult to apply this algorithm to the identification of an impact load or a random load.

Regularization and selection of the regularization parameter

Because of the inevitable ill-condition of the matrix $\mathit{H}$ , together with the noise contained in the measured response, the inverse problem is always ill-posed. In this article, Tikhonov regularization method is used to suppress the effect of noise on the identified results

$$minJ(\mathit{a},\lambda )=\Vert \mathit{Ha}-\mathit{y}{\Vert }^2+\lambda \Vert \mathit{a}{\Vert }^2$$

(17)

in which $\lambda$ is the regularization parameter and $\lambda >0$ , which is determined by the L-curve criterion.

After the regularization parameter is determined, the least-squares method is used to obtain the unknown sequence $\mathit{a}$

$$\mathit{a}=({\mathit{H}}^T\mathit{H}+\lambda \mathit{E}{)}^{-1}{\mathit{H}}^T\mathit{y}$$

(18)

in which $\mathit{E}$ is an identity matrix.

Next, the sequence $\mathit{a}$ calculated by equation (18) is substituted into equation (15) to obtain the load term and the parameter–load term, and then the unknown parameter perturbation terms can be obtained by the algorithm in the preceding section.

In the case of calculating the partial derivative matrix of unknown parameters in equation (8), it is necessary to take a very small perturbation step to ensure the accuracy of the final identified results, which also leads to the fact that the matrix ${\mathit{H}}_i$ $(i=1,2,\dots ,l)$ is much smaller than the matrix ${\mathit{H}}_0$ in the sense of norm. If determined directly using the L-curve criterion for equation (17), the regularization parameter is often very small, which makes the iterative results highly unsatisfactory.

Therefore, the perturbation terms in equation (17) are removed first

$$minJ(\mathit{a},\lambda )=\Vert {\mathit{H}}_0{\mathit{La}}_0-\mathit{y}{\Vert }^2+\lambda \Vert {\mathit{a}}_0{\Vert }^2$$

(19)

The regularization parameter is calculated using equation (19), which is then applied to equation (18). In addition, because the regularization parameter is calculated without considering ${\mathit{H}}_i$ $(i=1,2,\dots ,l)$ , its value is much larger than the result calculated directly using equation (17) and the unknown parameter perturbation terms calculated using this regularization parameter are naturally small, which leads to small changes in the identified values of the unknown parameters in each iteration step. Accordingly, the changes in ${\mathit{H}}_0$ and the regularization parameter between two adjacent iteration steps are also small, so there is no need to calculate the regularization parameter in each iteration step, thereby greatly improving the efficiency of the algorithm.

Procedures of the algorithm

The procedures for the simultaneous identification algorithm proposed in this article based on the perturbation method are as follows:

• Step 1. Measure the response of the structure with unknown parameters under the unknown load.

• Step 2. Assume initial values for the unknown parameters and calculate the impulse response functions between the response points and the load point and the first-order partial derivative of the impulse response functions with respect to all unknown parameters by the difference method. Then, form the impulse response matrix and partial derivative matrices.

• Step 3. Use the L-curve criterion to determine the regularization parameter, and then calculate the load term and the parameter–load term by Tikhonov regularization method.

• Step 4. Calculate the unknown parameter perturbation terms, and then after the determination of their positive and negative properties, add them to the initial values of the unknown parameters as the new initial values.

• Step 5. Repeat steps 2–4 until the identified results converge, but only after a certain number of iteration steps, recalculate the regularization parameter in step 3.

Numerical examples

A 10-storey shear frame

In this section, a 10-storey shear frame is chosen as the research model, as shown in Figure 1. $m_i=200\mathrm{kg}$ and $k_i=1000\mathrm{N}{\mathrm{m}}^{-1}$ $(i=1,2,\dots ,10)$ are the mass and stiffness of each storey, respectively. The structural damping is assumed to be proportional, and the mass matrix correlation coefficient $\alpha$ and the stiffness matrix correlation coefficient $\beta$ of the damping matrix are 0.04 s⁻¹ and 0.04 s, respectively. The load $f(t)$ acts on top of the frame, and the time history of the load is shown in Figure 2(f).

Figure 1.

A 10-storey shear frame.

Figure 2.

Identified results of one calculation in the case of 5% noise level and three response points: (a) convergence of identified $k_3$ , (b) convergence of identified $k_5$ , (c) convergence of identified $k_7$ , (d) convergence of identified $\alpha$ , (e) convergence of identified $\beta$ , and (f) comparison of identified $f(t)$ with the actual load.

In this article, the damage is simulated by reducing the stiffness value of the structure, and it is assumed that the structure has suffered 30% damage at $k_3$ , $k_5$ , and $k_7$ . Then, these three stiffness values are treated as unknown parameters. Because $\alpha$ and $\beta$ are difficult to determine accurately by simple tests, the two parameters are also treated as unknown parameters.

In this article, Newmark method ²⁸ is used to calculate the structural response. The time-sampling interval is 0.01 s, and the total calculation time is 10 s. Noise is added to the response of degrees of freedom chosen as the response points by the following equation

$${\mathit{y}}_n=\mathit{y}+L_n·\mathrm{std}(\mathit{y})·\mathrm{randn}$$

(20)

in which ${\mathit{y}}_n$ is the noise-polluted response, $L_n$ is the noise level, $\mathrm{std}(y)$ represents the standard deviation of the response, and randn is a Gaussian white noise sequence with a mean value of zero and a standard deviation of one.

Algorithm validation

First, the load and unknown parameters are identified without considering the effect of noise. The initial values of the stiffness are all assumed to be 1000 Nm⁻¹ and the initial values of $\alpha$ and $\beta$ are assumed to be 0.01 s⁻¹ and 0.01 s, respectively. The truncation ratio for this section is taken as 5%, and the response points take 1 and 3 degrees of freedom. The identified value and error of the load and the unknown parameters are shown in Table 1.

Table 1.

Identified results of the load and unknown parameters with no noise in the response.

RP	TR (%)	$k_3$		$k_5$		$k_7$		$\alpha$		$\beta$		f
		IV (Nm−1)	Error (%)	IV (Nm−1)	Error (%)	IV (Nm−1)	Error (%)	IV (s−1)	Error (%)	IV (s)	Error (%)	Error (%)
1, 3	5	700.00	0	700.00	0	700.00	0	0.0400	0	0.0400	0	3.144*e–5

RP: response points; TR: truncation ratio; IV: identified value.

As Table 1 shows, in the case of only two response points, the load and unknown parameters can be accurately identified, which verifies the correctness of the algorithm proposed in this article.

Effect of noise level and number of response points on the identified results

Although accurate identified results of the load and unknown parameters can be obtained without considering the effect of noise, the response always contains the noise components, and the existence of noise greatly amplifies the effect of the ill-conditioning of the inverse problem on the identified results. In this section, noise with different levels is added to the response to simulate the measured response, and different numbers of response points are taken to identify the load and unknown parameters to further verify the applicability of the proposed algorithm.

In addition, it should be noted that because the response contains random noise, even if the same level of noise is added to the response, the final identified results are certainly not the same, especially for unknown parameters. In this article, noise with the same level but with different distributions is added to the response, and multiple sets of identified results of the load and unknown parameters are obtained by multiple calculations. Finally, the mean values and standard deviations of the identified values of the unknown parameters are calculated. Because the mean value of an identified parameter is not necessarily equal to its actual value, the error of the ith unknown parameter is defined by the following equation

$$\mathrm{er}{r}_{b_i}=(|b_i^{\mathrm{actual}}-b_i^{\mathrm{mean}}|+b_i^{\mathrm{std}})/b_i^{\mathrm{actual}}(i=1,2,\dots ,l)$$

(21)

in which $b_i^{\mathrm{actual}}$ , $b_i^{\mathrm{mean}}$ , and $b_i^{\mathrm{std}}$ are the actual value, mean value, and standard deviation of the identified results of the ith unknown parameter $b_i$ , respectively.

The load identified error is determined by the average error of the multiple sets of the identified load. The results obtained in this way can better reflect the accuracy and stability of the algorithm.

The truncation ratio and initial values of the unknown parameters are the same as those in the preceding section. Table 2 shows the mean values and errors of the identified unknown parameters and the errors of the identified load in the case of different noise levels and different number of response points. The convergences of the unknown parameters and the comparison of the identified $f(t)$ and the actual load of one calculation in the case of 5% noise level and three response points (1, 3, and 5 degrees of freedom) are shown in Figure 2.

Table 2.

Identified results of the load and unknown parameters with different levels of noise and different number of response points.

RP	NL (%)	$k_3$		$k_5$		$k_7$		$\alpha$		$\beta$		$f$
		MV (Nm−1)	Error (%)	MV (Nm−1)	Error (%)	MV (Nm−1)	Error (%)	MV (s−1)	Error (%)	MV (s)	Error (%)	Error (%)
1, 3	1	700.18	0.16	700.18	0.18	700.27	0.21	0.0401	1.24	0.0399	0.76	1.56
	5	700.73	0.68	700.67	0.87	700.42	1.08	0.0402	5.13	0.0397	2.97	5.16
	10	700.55	1.18	701.19	1.37	700.78	1.68	0.0403	7.66	0.0396	5.59	7.56
1, 3, 5	1	699.93	0.12	700.09	0.14	700.14	0.16	0.0400	1.06	0.0399	0.67	1.42
	5	700.35	0.43	700.21	0.48	700.02	0.55	0.0402	3.18	0.0398	2.31	5.02
	10	700.51	0.84	700.53	0.89	699.74	0.94	0.0401	5.54	0.0397	3.85	7.32
1, 3, 5, 7	1	700.01	0.10	700.14	0.13	700.09	0.13	0.0400	0.90	0.0399	0.61	1.40
	5	700.19	0.36	700.37	0.40	700.17	0.43	0.0401	2.72	0.0398	1.84	4.91
	10	699.56	0.73	700.27	0.77	700.31	0.78	0.0401	4.34	0.0396	3.36	7.16

RP: response points; NL: noise level; MV: mean value.

In summary, satisfactory identified results of the load and unknown parameters are obtained in all cases in Table 2, especially the structural stiffness. Even in the case of 10% noise levels, the errors of the identified stiffness are not greater than 1.7% when only two response points are used, which is completely acceptable in engineering.

Furthermore, as the noise level increases, the errors of the identified results increase, and by increasing the number of response points, the errors decrease significantly. Therefore, when more accurate identified results are necessary, more response points must be selected, especially when the response is seriously polluted by noise. In addition, the convergences of unknown parameters have obvious convergence points as shown in Figure 2(a)–(e); thus, there is no requirement for a convergence condition for the proposed algorithm.

Effect of truncation ratio on the identified results

A parameter called the truncation ratio is defined in this article; thus, this section studies its effect on the identified results.

To compare the effect of the truncation ratio on the identified results more directly and effectively, in the case of 5% noise level and three response points (1, 3, and 5 degrees of freedom), multiple sets of responses with the same noise level, but the noise is different in the time history, are generated first. Next, the multiple generated sets of responses are used to identify the load and unknown parameters under the different truncation ratios from 0% to 10%. Finally, the mean values and errors of the identified results are obtained at different truncation ratios. The initial values of the unknown parameters are the same as those in the preceding section. The curves of mean values and errors of the identified results of the stiffness in the case of different truncation ratios are shown in Figure 3.

Figure 3.

Identified results of stiffness with different truncation ratios: (a) mean value and (b) error.

As Figure 3 shows, when the truncation ratio is 0%, that is, when the positive and negative properties of the perturbation terms of the unknown parameter are determined in the case where the sequence $\mathit{r}{\mathit{f}}_i$ is not truncated, the errors of the identified stiffness are very large, and the mean values also have large errors compared to the actual value. When the truncation ratio increases to 1%, the errors obviously decrease, and the mean values become close to the actual value. In the range of 2%–8% of the truncation ratio, the errors are very small and stable. When the truncation ratio is greater than 8%, the error of the identified results will increase.

Therefore, a conclusion can be drawn that the introduction of the truncation ratio effectively reduces the error and increases the stability of the proposed algorithm. In addition, in a very large range of the truncation ratio, the identified value can approach the actual value, illustrating that the truncation ratio is not a difficult parameter to determine.

Effect of initial values and norm types on the identified results

To study the effect of initial values of unknown parameters on the identified results, with different initial values of the stiffness, the convergences of $k_3$ are shown in Figure 4(a). In addition, in view of the norm type mentioned in calculating the unknown parameter perturbation terms, the convergences of $k_3$ with different norm types are shown in Figure 4(b).

Figure 4.

Convergences of identified $k_3$ with (a) different initial values and (b) different norm types.

As Figure 4 shows, with different initial values and different norm types, the final identified values of $k_3$ are basically the same, illustrating that the proposed algorithm is not sensitive to the above two factors.

A simply supported beam

A simply supported beam with a uniform cross-section is used as the research model in this section, as shown in Figure 5. The transverse vibration equation of a beam under a concentrated load is as follows

$$\begin{array}{c}\rho A\frac{{\partial }^{2}w(x,t)}{\partial t^2}+c\frac{\partial w(x,t)}{\partial t}+\mathrm{EI}\frac{\partial w(x,t)}{\partial x^4}\hfill \\ =f(x,t)=\delta (x-x_p)p(t)\hfill \end{array}$$

(22)

in which $\rho$ is the density, A is the area of the cross-section, c is the damping, E is the elastic modulus, I is the inertia moment of the cross-section with respect to the neutral axis, $\omega (x,t)$ is the transverse vibration displacement, $f(x,t)$ is the external load, $p(t)$ is the time history of the load, and $x_p$ is the load position.

Figure 5.

A simply supported beam.

With the basic idea of the reference, ²⁹ $p(t)$ is expanded with respect to the Legendre polynomials

$$\begin{array}{c}\multicolumn{1}{c}{p(t)=\sum _{v=1}^V{a}_v{L}_v(t)=\left[\begin{array}{cccc}{L}_1\hfill & L_2\hfill & \cdots \hfill & L_V\hfill \end{array}\right]}\\ {\left[\begin{array}{cccc}{a}_1\hfill & a_2\hfill & \cdots \hfill & a_V\hfill \end{array}\right]}^T\hfill \end{array}$$

(23)

Using the orthogonality of the mass, damping, and stiffness with respect to the structural modal shapes, one can obtain

$$\{\begin{array}{c}{M}_i=\underset{0}{\overset{L}{\int }}\rho A{{\varphi }_i}^2(x)\mathrm{dx}\hfill \\ \multicolumn{1}{c}{C_i=\underset{0}{\overset{L}{\int }}c{{\varphi }_i}^2(x)\mathrm{dx}}\\ \multicolumn{1}{c}{K_i=\underset{0}{\overset{L}{\int }}\mathrm{EI}{({{\varphi }_i}^{\prime \prime }(x))}^2\mathrm{dx}(i=1,2,\dots ,\infty )}\\ \multicolumn{1}{c}{P_i(t)=\underset{0}{\overset{L}{\int }}\delta (x-x_p)p(t){\varphi }_i(x)\mathrm{dx}=p(t){\varphi }_i(x_p)}\end{array}$$

(24)

in which ${\phi }_i(x)$ , $M_i$ , $C_i$ , $K_i$ , and $P_i(t)$ are the ith order modal shape, modal mass, modal damping, modal stiffness, and modal load, respectively, and L is the length of the beam.

Substituting equation (23) into the expression of the modal load in equation (24), one can obtain

$$P_i(t)={\varphi }_i(x_p)[C_1(t),C_2(t),\cdots ,C_V(t)][a_1,a_2,\cdots a_V{]}^T$$

(25)

The dynamic response of the simply supported beam can also be obtained using Duhamel integral together with the modal superposition method

$$\omega (x,t)=\sum _{i=1}^{\infty }\left(\frac{{\varphi }_i(x){\varphi }_i(x_p)}{M_i{\omega }_{\mathrm{di}}}{\left\{\begin{array}{c}\underset{0}{\overset{t}{\int }}{h}_i(t-\tau )L_1(\tau )d\tau \hfill \\ \multicolumn{1}{c}{\underset{0}{\overset{t}{\int }}{h}_i(t-\tau )L_2(\tau )d\tau }\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\underset{0}{\overset{t}{\int }}{h}_i(t-\tau )L_V(\tau )d\tau }\end{array}\right\}}^T\right)\left\{\begin{array}{c}{a}_1\hfill \\ \multicolumn{1}{c}{a_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{a_V}\end{array}\right\}$$

(26)

in which $h_i(t)=e^{-{\xi }_i{\omega }_{i}t}\sin ({\omega }_{\mathrm{di}}t)$ is the ith order impulse response function, ${\omega }_{\mathrm{di}}={\omega }_i\sqrt{1-{{\xi }_i}^2}$ is the ith order damped natural frequency, ${\xi }_i=C_i/(2M_i{\omega }_i)$ is the ith order modal damping ratio, and ${\omega }_i=\sqrt{K_i/M_i}$ is the ith order natural frequency.

Suppose that $N_d$ response points, Q time-sampling points, and the previous $N_m$ orders of the modal shapes are taken. The response can be obtained as follows

$$\begin{array}{c}\left\{\begin{array}{c}\omega (x_1,t_1)\hfill \\ \multicolumn{1}{c}{\omega (x_1,t_2)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\omega (x_1,t_Q)}\\ \multicolumn{1}{c}{\omega (x_2,t_1)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\omega (x_{N_d},t_Q)}\end{array}\right\}=\hfill \\ \left(\sum _{i=1}^{N_m}\frac{{\varphi }_i(x_p)}{M_i{\omega }_{\mathrm{di}}}\left\{\begin{array}{c}{\varphi }_i(x_1)\left[\begin{array}{cccc}{H}_{i,1}^{t_1}\hfill & H_{i,2}^{t_1}\hfill & \cdots \hfill & H_{i,V}^{t_1}\hfill \end{array}\right]\hfill \\ \multicolumn{1}{c}{{\varphi }_i(x_1)\left[\begin{array}{cccc}{H}_{i,1}^{t_2}\hfill & H_{i,2}^{t_2}\hfill & \cdots \hfill & H_{i,V}^{t_2}\hfill \end{array}\right]}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\varphi }_i(x_1)\left[\begin{array}{cccc}{H}_{i,1}^{t_Q}\hfill & H_{i,2}^{t_Q}\hfill & \cdots \hfill & H_{i,V}^{t_Q}\hfill \end{array}\right]}\\ \multicolumn{1}{c}{{\varphi }_i(x_2)\left[\begin{array}{cccc}{H}_{i,1}^{t_1}\hfill & H_{i,2}^{t_1}\hfill & \cdots \hfill & H_{i,V}^{t_1}\hfill \end{array}\right]}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\varphi }_i(x_{N_d})\left[\begin{array}{cccc}{H}_{i,1}^{t_Q}\hfill & H_{i,2}^{t_Q}\hfill & \cdots \hfill & H_{i,V}^{t_Q}\hfill \end{array}\right]}\end{array}\right\}\right)\left\{\begin{array}{c}{a}_1\hfill \\ \multicolumn{1}{c}{a_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{a_V}\end{array}\right\}\hfill \end{array}$$

(27)

in which $H_{i,v}^{t_q}={\int }_0^t{h}_i(t_q-\tau )L_v(\tau )d\tau$ , $i=1,2,\dots ,N_m$ , $q=1,2,\dots ,Q$ , and $v=1,2,\dots ,V$ .

The above equation can be simplified as

$${\left\{\omega \right\}}_{(N_d\times Q)\times 1}=[H{]}_{(N_d\times Q)\times V}·{\left\{a\right\}}_{V\times 1}$$

(28)

Thus, the same form as equation (2) is obtained, and the load and structural unknown parameters can be identified simultaneously by the algorithm proposed in the preceding section.

In this article, the cross-section of the beam is rectangular; the beam height, width, and length are 0.005, 0.06, and 2 m, respectively; the elastic modulus E is 1.8e11 Nm⁻²; the density $\rho$ is 7850 kgm⁻³; and the damping c is 3 Nsm⁻¹. The load acts at a position of 0.7 m along the beam, and the time history of $p(t)$ is shown in Figure 6(d).

Figure 6.

Identified results of one calculation in the case of set 4: (a) convergence of identified E, (b) convergence of identified c, (c) convergence of identified $x_p$ , and (d) comparison of identified $p(t)$ with the actual load.

The elastic modulus E, damping c, and load position $x_p$ are considered to be unknown. The initial values of the elastic modulus E and the damping c are assumed to be 2e11 Nm⁻² and 1 Nsm⁻¹, respectively; the initial position of the load $x_p$ is assumed to be at 1 m along the beam. The time-sampling interval is 0.001 s, and the total calculation time is 1 s. The noise level and the truncation ratio are both set to 5%.

To study the effect of the distribution of the response points on the identified results, six sets of response points are selected, and the identified results of the load and unknown parameters with different sets of response points are shown in Table 3. The convergences of the unknown parameters and the comparison of the identified $p(t)$ and the actual load of one calculation in the case of set 4 are shown in Figure 6.

Table 3.

Identified results of the load and unknown parameters with different sets of response points.

Set no.	RP (m)	E		c		$x_p$		p
		MV (Nm−2)	Error (%)	MV (Nsm−1)	Error (%)	MV (m)	Error (%)	Error (%)
1	0.85, 0.95, 1.05, 1.15	1.795e11	1.18	2.967	5.04	0.7024	1.21	5.40
2	0.70, 0.90, 1.10, 1.30	1.798e11	0.55	2.976	3.79	0.7011	0.46	4.56
3	0.65, 0.85, 1.15, 1.35	1.797e11	0.54	2.993	2.56	0.7018	0.43	4.54
4	0.40, 0.80, 1.20, 1.60	1.796e11	0.55	2.997	2.11	0.7013	0.33	4.52
5	0.25, 0.75, 1.25, 1.75	1.796e11	0.57	2.991	2.78	0.7013	0.35	4.81
6	0.10, 0.70, 1.30, 1.90	1.803e11	0.67	3.011	3.37	0.6989	0.49	4.94

RP: response points; MV: mean value.

As Table 3 shows, the more concentrated the distribution of response points, the larger the errors, even though the response points are distributed in the vicinity of the midpoint of the beam, where the response is relatively large. This is because the response history of one point is too much like the others under these circumstances, which leads to a reduction in the information about the load and structural parameters contained in the response. When two of the four response points approach the two endpoints of the beam, the errors are slightly larger. Therefore, it is also not recommended for the response points to be distributed near the two endpoints of the beam because the response of this part of the beam is much smaller.

Conclusion

This article proposes a new algorithm for the simultaneous identification of the load and unknown parameters of the structure using only a few response points based on the perturbation method. The algorithm only requires initial values for the unknown parameters, which are not required for the loads. In addition, the proposed algorithm does not require a convergence condition for the iteration process. In addition, the regularization parameter does not need to be calculated in each iteration step, thereby effectively increasing the efficiency of the algorithm. Furthermore, a parameter called the truncation ratio is defined, and the identified results for an example show that the introduction of this parameter greatly reduces the errors obtained with the proposed algorithm and that the best range of this parameter is relatively large. In a case with different noise levels and different numbers of response points, satisfactory results are all obtained, illustrating that the algorithm provides good accuracy and applicability. Moreover, the proposed algorithm is also applied to the identification of a simply supported beam, and the study of the distribution of response points shows that when the response points are distributed over the full beam, the identified results are better than that when the response points are distributed over only part of the beam.