# Identification of Sudden Stiffness Change in the Acceleration Response of a Nonlinear Hysteretic Structure

Authors: Sheng-Lan Ma; Shao-Fei Jiang; Chen Wu; Si-Yao Wu
Journal: Shock and Vibration
 (2020)

---

## Abstract

The integration of discrete wavelet transform and independent component analysis (DWT-ICA) method can directly identify time-varying changes in linear structures. However, better metrics of structural seismic damage and future performance after an event are related to structural permanent and total plastic deformations. This study proposes a two-stage technique based on DWT-FastICA and improved multiparticle swarm coevolution optimization (IMPSCO) using a baseline nonlinear Bouc–Wen structural model to directly identify changes in stiffness caused by damage as well as plastic or permanent deflections. In the first stage, the measured structural dynamic responses are preprocessed firstly by DWT, and then the Fast ICA is used to extract the feature components that contain the damage information for the purpose of initially locating damage. In the second stage, the structural responses are divided at the identified damage instant into segments that are used to identify the time-varying physical parameters by using the IMPSCO, and the location and extent of damage can accordingly be identified accurately. The efficiency of the proposed method in identifying stiffness changes is assessed under different ground motions using a suite of two different ground acceleration records. Meanwhile, the effect of noise level and damage extent on the proposed method is also analyzed. The results show that in a realistic scenario with fixed filter tuning parameters, the proposed approach identifies stiffness changes within 1.25% of true stiffness within 8.96 s; therefore, it can work in real time. Parameters are identified within 14% of the actual as-modeled value using noisy simulation-derived structural responses. This indicates that, in accordance with different demands, the proposed method can not only locate and quantify damage within a short time with a high precision but also has excellent noise tolerance, robustness, and practicality.


---

## Body

## 1. Introduction

Civil engineering structures are subjected to continuous structural deterioration caused by aging, low-cycle fatigue loads from smaller earthquakes, and daily environmental loading. Therefore, the severity of deterioration needs to be monitored periodically in order to ensure structural integrity and safety [1]. If damage is detected in the early stage, maintenance works will be carried out timely with low cost. For this purpose, structural responses (e.g., accelerations) due to structural vibrations induced by external loading are on-time collected by structural health monitoring (SHM) systems; analyzing and processing collected structural response is a popular way to identify, locate, and quantify structural damage based on the principle that damage affects the mechanical properties of the structure (i.e., stiffness and damping), which will change structural dynamic properties (e.g., frequencies and model shapes) [2, 3].A further approach to damage detection is the identification of an anomaly in the time history. A number of researchers have been able to locate an anomaly in the processed signal that is related to damage by using signal processing techniques [4], which often combine frequencies and time domains (such as wavelet analysis [5], empirical mode decomposition, and blind source separation). These approaches mainly aim to detect the anomaly in online or real-time structural health monitoring. For example, Yang and Nagarajaiah [6] proposed a unsupervised blind source separation technique, which combined wavelet transform and independent component analysis (ICA), for the sake of detecting damage instant for linear structure. In addition, Jiang et al. [7] employed wavelet transform combined with data fusion to identify structural damage anomaly in time domain for linear structures. Generally, the majority of proposed damage instant detection methods are used to identify linear structures. Unfortunately, when structures are subjected to earthquake excitation, nonlinearity widely exists in these damaged structures [8]. For instance, when damage occurs in reinforced concrete structures, cracks opening and closing under dynamic excitation are typical nonlinear process which leads to the hysteretic performance of the structures [9].After completing the structural response novelty detection in time history due to sudden stiffness change, it is necessary to use an accurate and fast approach to identify the degree of the stiffness changed and the parameters in nonlinear hysteretic system for urgent postevent response and decision making. Therein, the particle swarm optimization (PSO) algorithm, due to its excellent characteristics of simple principle and fast convergence, is widely applied in system identification [10], damage detection [11, 12], and sensor optimization placement [13] in linear structure, as well as in parameters identification [14] in nonlinear hysteretic structure. The principle of optimization-based identification consists in minimizing the difference between experimentally and numerically obtained structural dynamic properties. However, the standard PSO algorithm may get trapped in the local optimum and might not converge to the global optimum. In order to improve the performance of standard PSO algorithm, some investigation was made by combining PSO algorithm with other algorithm and modifying the inertia weights. By introducing some mechanism into standard PSO algorithm, Charalampakis and Dimou [15] presented two variants of PSO algorithm to identify the Bouc–Wen hysteretic system. Furthermore, Zhang and Xia [16] developed an improved PSO algorithm by defining a fitness function to identify parameters of nonlinear dynamic hysteretic models. Overall, the fitness functions for parameters identification based on PSO algorithm mainly utilized modal information. However, due to the significant computational complexity and longer runtime, the applications of aforementioned methods for real-time or time varying damage detection are limited.Due to the fact that few methods can not only detect structural response novelty but also identify structural damage severity in nonlinear structures, this paper focuses on developing an approach that will capture a nonlinear behavior in time domain to detect structural response novelty firstly and then utilize an optimization algorithm, with fitness function in time domain, to identify structural parameters and structural stiffness in less time. The first question that arises is how to characterize a nonlinear response as mentioned above. This has been addressed in a wide range of mechanical and civil engineering applications via signal processing methods like Hilbert–Huang transform-based [1, 17], Wavelet transform-based [7], neural network-based [18, 19], and independent component analysis-based methods [20]. Yang and Nagarajaiah [6] combined wavelet transform with independent component analysis (ICA) to detect linear structural instant, which had excellent robustness. As the extension of Reference [6], this paper firstly utilized the integration of discrete wavelet transform and independent component analysis (DWT-ICA) to detect nonlinear structural response novelty. The second question that arises is how to quickly and accurately identify nonlinear structural stiffness and finally to achieve the aim of detecting sudden stiffness change. Just as the aforementioned study and the authors’ previous study [21], an improved multiparticle swarm coevolution optimization algorithm (IMPSCO) with time-domain fitness function is applied to detect nonlinear structural damage by the first time. All in all, this paper presents a two-stage technique based on DWT-FastICA and the improved multiparticle swarm coevolution optimization (IMPSCO) using a baseline nonlinear Bouc–Wen structural model to directly identify sudden stiffness changes in the acceleration response.The organization of the paper is as follows. Section2 reviews the basic theory of the Bouc–Wen model. Section 3 numerically and experimentally demonstrates the feasibility of IMPSCO algorithm to identify time-varying nonlinear structural parameters. Section 4 depicts the damage detection strategy based on DWT-FastICA and IMPSCO algorithm. Numerical simulations of two-DOF nonlinear hysteretic system are carried out in Section 5, and Section 6 gives the concluding remarks finally.


## 2. Bouc–Wen Model: Model Formulation and Parameter Constraints

Structural dynamic behavior must generally be taken into account in the design of mechanical systems to insure their performance and reliability [22]. In structural systems, hysteresis appears as a natural mechanism in the materials used and produces restoring forces that dissipate energy. Hysteresis, as used here, refers to the memory nature of inelastic structural behaviors where the restoring force depends not only on the instantaneous deformations but also on the history of the deformations. The detailed analytical modeling of this behavior results in very complicated nonlinear models that are not suitable for sudden stiffness change. One of the most common semiphysical models proposed is the first-order nonlinear differential equation known as the Bouc–Wen model. More specifically, a nonlinear frame structure with passive MR dampers under seismic base motion, as shown in Figure 1, can be modeled as follows [23]:(1)MX¨t+CX˙t+RX,Z,t=−Ft,where M and C are the mass and damping matrices of the system, respectively; X, X˙, and X¨ are displacement, velocity, and acceleration matrix of the base-isolation system, respectively; F is the external force; and R is the vector of restoring force governed by the following first-order differential equation from the so-called classical Bouc–Wen model:(2)RX,Z,t=R1−R2,R2−R3,⋯,RN−1−RNT,in which(3)Ri=αikiXit+1−αikiYiZit,Z˙it=1YiAiX˙it−βiX˙itZitZitni−1+γiX˙itZitni,Ai>0,βi>0,−βi<γi<βi,ni≥1,where k is the stiffness of the system; A, β, γ, and n are stiffness, loop fatness, loop pinching, and abruptness parameters in the classical Bouc–Wen model, respectively; α is the bilinear factor, defined as the ratio of the post- to preyield stiffness of the system; Y is the yield displacement of the system; Zt is the dimensionless Bouc–Wen hysteresis component; the subscript i represents the i-th story; and N is the total number of stories.Figure 1
Model of a base-isolation system with passive MR dampers.It is noted thatn is the power factor which determines the sharpness of the curve from elastic to plastic force-deflection behavior of the system.In this paper, we assume structural state vectorS as follows:(4)S=X,X¨,Z.Therefore, equation (1) can be transformed as a first-order ordinary differential equation as follows:(5)S=X,X¨,Z=fS.Therein,(6)X¨=M−1Ft−CX˙t−RX,Z,t.As a result, if we initialized structural state vector, equation (6) can be solved by the fourth-order Runge–Kutta method [24] to calculate structural dynamic responses.


## 3. Feasibility of IMPSCO Algorithm Applied in Identification of Nonlinear Structural Parameters

In order to identify time-varying nonlinear structural parameters, the improved multiparticle swarm coevolution optimization (IMPSCO) algorithm is applied, which is further used to locate and quantify structural damage as well as identify the parameters in the Bouc–Wen model.


### 3.1. Theory of the IMPSCO Algorithm

For basic multiparticle swarm coevolution optimization algorithm (MPSCO) [25], multiple subpopulations are divided into two layers. All particles from the upper layer follow the optimum of the entire population so as to obtain a faster convergence speed, while all particles from the lower layer follow the optimum of the subpopulation to ensure the population diversity. Although the performance of basic MPSCO is better than standard PSO in some respects, the subpopulations in the lower layer still perform the process of standard PSO, which makes particle falling into the local optimum possible. To solve this problem, we proposed the IMPSCO algorithm, and the details can be seen in Reference [21], and it was applied to locate and quantify damage in a linear structure. Figure 2 shows the flowchart of the IMPSCO algorithm.Figure 2
The flowchart of the IMPSCO algorithm.Wherein, the key point in IMPSCO is worst particle replacement. More specially, when the particle is recorded as the worst for the predetermined timesIw, replace it with local center of gravity Gg:(7)Gg=∑i=1szis,where s and zi represent the number of the selected excellent particles and their position, respectively.The most important task for optimization is to determine the fitness function in the IMPSCO. A useful fitness function not only has high accuracy but also needs less runtime, which are the main evaluation indices of the optimization algorithm for detecting sudden structural stiffness change. In addition, some parameters in the nonlinear hysteretic model may have little effects on the shape of the hysteretic curve. Therefore, we utilized the difference between the calculated and testing accelerations as the fitness function value, which contains lots of structural sensitive information. However, it is difficult to calculate the minimum. We employed the reciprocal of the mean square error (MSE) between the reproduced acceleration responsesX¨tp and the reference measured acceleration responses X¨t and replaced the mean square error (MSE) between them as the fitness function. When cast in discrete form, it can be expressed as(8)maxfθ=1∑i=1N∑kl=1LX¨itkl−X¨itklp2,where θ is the structural parameters vector to be optimized; N is the total number of measuring points; and L is the length of time history.


### 3.2. Feasibility: Numerical Simulation

#### 3.2.1. Model Description

To evaluate the performance of IMPSCO algorithm applied to identification of nonlinear hysteretic base-isolation systems, a realistic base-isolated system is created as shown in Figure3 [26]. The numerical model can be simplified as a two degree-of-freedom (DOF) structure. Its basic definition includes the following: mass of each floor: m1=m2=1kg; initial stiffness of each floor: k1=k2=102.64N/m; damping coefficients of each floor: ζ1=ζ2=2%; parameters in the classical Bouc–Wen model: A=1;β=0.5;γ=0.5;n=2;α=0.1; the yield displacement of the system: Y=0.045m; and the basic period of structure: T=1s.Figure 3
2-DOF structure with Bouc–Wen model subjected to seismic excitation.Meanwhile, according to the sensitivity analysis of the parameters in the classical Bouc–Wen model listed in Reference [26], we assumed that some parameters with low sensitivity are constants in some cases in order to reduce the identified parameters in IMPSCO and simply the identification process. In summary, this study considered four cases with different parameters known. In addition, the search range of each parameter (θ) is within the interval of [0.5b, 2b], where b is the theoretical value of the corresponding parameter (θ).In this study, the structure was subjected to El-Centro earthquake wave, Kobe earthquake wave, and Northridge earthquake wave, respectively, at the base, as shown in Figure3. In order to reduce the particle dimension of IMPSCO, four cases were simulated to observe the feasibility of IMPSCO algorithm, as shown in Table 1. Therefore, by setting the sampling frequency of 50 Hz and the loading time t = 30 s, the structural acceleration responses were calculated by using the classical fourth-order Runge–Kutta algorithm [24].Table 1
Different cases for simulation.
Cases
Unknown parameters
Known parameters
S1
k
1
,
k
2
,
ζ
1
,
ζ
2
,
Y
,
α
,
β
,
γ
,
n
—
S2
k
1
,
k
2
,
ζ
1
,
ζ
2
,
Y
,
α
,
β
,
n
γ
=
0.5
S3
k
1
,
k
2
,
ζ
1
,
ζ
2
,
Y
,
β
,
n
α
=
0.1
;
γ
=
0.5
S4
k
1
,
k
2
,
ζ
1
,
ζ
2
,
Y
,
α
,
β
,
γ
n
=
2To consider the operating environment, Gaussian white noise with different levels was added to all structural responses and the earthquake motion. The structural responses were then represented as(9)yi=yia×1+εR,where yi and yia represent the contaminated and theoretical signals, respectively; R is a normally distributed random variable with zero mean and a derivation of 1; and ε is an index representing the noise level.Therefore, the signal-to-noise ratio (SNR) is defined as(10)SNR=20log101ε=20log10AsignalAnoise,where Asignal and Anoise are the amplitudes of the signal and noise, respectively.In this paper, the SNR was 30 dB.


#### 3.2.2. Identification Process

In order to identify parameters of 2-DOF structure, two issues need to be ensured for the IMPSCO algorithm: encoding parameters and initial setting parameters. The details can be seen as follows.(1) Parameter Encoding Involved in the IMPSCO. This study assumed that the mass of each floor was known. Therefore, the parameters encoded in the IMPSCO were initial stiffness (k1, k2) and damping coefficients (ζ1, ζ2), yield displacement (Y) of the first floor, bilinear factor (α), loop fatness (β), loop pinching (γ), and abruptness parameters (n) in the classical Bouc–Wen model system, and the details could be seen in Figure 4. Consequently, the number of subpopulations m was 9.Figure 4
The parameters encoded in the IMPSCO.(2) Parameters in IMPSCO. To maximize the fitness function using IMPSCO presented in Section 3.1, some parameters should be set and initialized in advance for IMPSCO. Generally, 3∼5 subpopulations with 30 to 100 particles in total were sufficient. The inertia weight w was set linearly varying from 0.9 to 0.4 and learning factors c1 and c2 were set to 2 simultaneously. In addition, the limited times for the worst record Iw and the number of the selected excellent particles s were both empirically set to 5∼10, respectively. Finally, the parameters of IMPSCO in this study were initialized as shown in Table 2.Table 2
Specifications of IMPSCO algorithm.
Parameter
Description
Value
mm
The total number of subpopulations
3
nn
The size of a subpopulation
10
c
1
Cognitive parameter
2
c
2
Social parameter
2
ρ
min
Minimum of inertia weight
0.4
ρ
max
Maximum of inertia weight
0.9
I
w
Limited times for worst record
5
S
The number of the selected excellent particles
6
I
max
The maximum iteration times
200


#### 3.2.3. Identification Results and Discussion

According to the process described in Section3.1, the responses shown in Figure 5 were employed to identify the parameters in the nonlinear structure and in Bouc–Wen model system. The results are listed in Table 3, and relative errors between identified results and theoretical values are shown in Table 4.Figure 5
Acceleration responses of the simulated structure under different earthquake waves,∗a1 is the acceleration responses of the first floor; a2 is the acceleration responses of the second floor; aq is corresponding earthquake excitation. (a) El-Centro earthquake waves. (b) Kobe earthquake waves. (c) Northridge earthquake waves.
(a)
(b)
(c)Table 3
Identification results of the simulated structure subjected to the different earthquake excitation in different cases.
Case
Seismic excitation
k
1 (N/m)
k
2 (N/m)
ζ
1 (%)
ζ
2 (%)
Y (mm)
α
β
γ
n
S1
(a) El-Centro earthquake
102.45
102.73
1.86
2.06
46.46
0.09
0.55
0.52
2.06
(b) Kobe earthquake
102.79
102.73
2.12
1.83
47.38
0.10
0.56
0.57
1.95
(c) Northridge earthquake
103.01
102.58
1.96
1.98
49.60
0.10
0.61
0.61
1.97
S2
(a) El-Centro earthquake
103.15
102.46
2.09
2.16
44.13
0.11
0.47
—
2.01
(b) Kobe earthquake
102.65
102.74
2.12
1.82
44.65
0.10
0.49
—
1.96
(c) Northridge earthquake
101.46
103.09
2.38
2.09
45.19
0.10
0.46
—
2.07
S3
(a) El-Centro earthquake
102.86
102.70
2.03
2.24
44.70
—
0.48
—
2.01
(b) Kobe earthquake
103.13
102.71
2.15
1.84
44.39
—
0.48
—
1.93
(c) Northridge earthquake
102.72
103.01
2.28
2.05
44.71
—
0.46
—
1.97
S4
(a) El-Centro earthquake
103.06
102.55
2.08
2.17
47.78
0.11
0.55
0.61
—
(b) Kobe earthquake
102.54
102.73
2.16
1.84
47.62
0.10
0.56
0.58
—
(c) Northridge earthquake
102.94
102.65
2.16
2.01
46.23
0.10
0.52
0.56
—Table 4
The relative errors between identified parameters and actual model parameters.
Case
Seismic excitation
Relative error (%)
k
1
k
2
ζ
1
ζ
2
Y
α
β
γ
n
S1
(a) El-Centro earthquake
0.19
0.09
7.00
3.00
3.24
10.00
10.00
4.00
3.00
(b) Kobe earthquake
0.15
0.09
6.00
8.50
5.29
0.00
12.00
14.00
2.50
(c) Northridge earthquake
0.36
0.06
2.00
1.00
10.22
0.00
22.00
22.00
1.50
Average values
0.23
0.08
5.00
4.17
6.25
3.33
14.67
13.33
2.33
S2
(a) El-Centro earthquake
0.50
0.18
4.50
8.00
1.93
10.00
6.00
—
0.50
(b) Kobe earthquake
0.01
0.10
6.00
9.00
0.78
0.00
2.00
—
2.00
(c) Northridge earthquake
1.15
0.44
19.00
4.50
0.42
0.00
8.00
—
3.50
Average values
0.55
0.24
9.83
7.17
1.04
3.33
5.33
—
2.00
S3
(a) El-Centro earthquake
0.21
0.06
1.50
12.00
0.67
—
4.00
—
0.50
(b) Kobe earthquake
0.48
0.07
7.50
8.00
1.36
—
4.00
—
3.50
(c) Northridge earthquake
0.08
0.36
14.00
2.50
0.64
—
8.00
—
1.50
Average values
0.26
0.16
7.67
7.50
0.89
—
5.33
—
1.83
S4
(a) El-Centro earthquake
0.41
0.09
4.00
8.50
6.18
10.00
10.00
22.00
—
(b) Kobe earthquake
0.10
0.09
8.00
8.00
5.82
0.00
12.00
16.00
—
(c) Northridge earthquake
0.29
0.01
8.00
0.50
2.73
0.00
4.00
12.00
—
Average values
0.27
0.06
6.67
5.67
4.91
3.33
8.67
16.67
—The following conclusions and remarks can be made from the two tables:(1)
The relative errors of initial stiffness of system (k) and the abruptness parameter (n) in the classical Bouc–Wen model are less than 3.5%, and the rest of identified parameters have lower computational accuracy. More specifically, the maximum relative errors for stiffness (α), loop fatness (β), and loop pinching (γ) are 10.00% (S1, S2, S4), 22.00% (S1), and 22.00% (S1, S4), respectively; the corresponding values are 10.22% (S1), 1.93% (S2), 1.36% (S3), and 6.18% (S4) for the yield displacement (Y); in addition, the maximum relative errors for damping coefficients (ζ1 and ζ2) are 8.50% (S1), 19.00% (S2), 14.00% (S3), and 8.50% (S4), respectively.There are two reasons for this phenomenon. The first one is that when we used the error between reproduced and measured acceleration responses as the fitness function, the stiffness of the system mainly influences the reproduced accelerations. By contrast, the damping coefficients have little effect on the reproduced accelerations. As a result, it is difficult to accurately identify the damping coefficients. The second reason is that, as presented in Reference [27], the abruptness parameter in the classical Bouc–Wen model (n) is a power factor, which changes the smoothness of transition from elastic region to plastic region in the hysteretic versus actual displacement curves, with little change in the acceleration response. Correspondingly, the remaining three parameters in the classical Bouc–Wen model are not sensitive to accelerations.(2)
It is accurate to identify stiffness of the system (k1,k2) and the abruptness parameter (n) in the classical Bouc–Wen model by using IMPSCO with the parameter listed in Table 2. Furthermore, for four cases (with different parameters to be identified and different particle dimension of IMPSCO), all the parameters of the nonlinear structure can be basically identified by IMPSCO. Although the identified results for the parameters in the classical Bouc–Wen model have lower accuracy than that of stiffness, the IMPSCO algorithm is suitable for identification of sudden stiffness change in actual structure. Just as shown in Table 4, the overall ranking for stiffness (k1,k2) and damping coefficients (ζ1,ζ2) is S1 > S4 > S3 > S2, and the overall ranking for remaining parameters is S3 > S2 > S4 > S1. This indicates that for all the parameters unknown case (S1), the identified results for the stiffness and damping have excellent accuracy. By contrast, the identified results for the parameters in the Bouc–Wen model have lower accuracy than other cases. However, the maximum error in the identification process is less than 16.7%. These error values are well within construction errors. It also demonstrates that the IMPSCO algorithm is applicable for detecting structural damage induced by stiffness change in structures.


### 3.3. Feasibility: Experiment

#### 3.3.1. Test Setup

A MR damper test [27] was presented to further validate the proposed IMPSCO algorithm. The test’s configuration schematic and the photograph of the testing system are shown in Figure 6.Figure 6
MR damper performance test setup [27].It is noted that the modified Bouc–Wen model is superior in describing the force-displacement relationship of MR damper. Its expression is as follows:(11)F=cx˙+kx+αz+f0,z˙=Ax˙−βx˙zzn−1+γx˙zn,where F and f0 are the damping force and initial damping force; c and k are damping and the stiffness of the system, respectively; α is the bilinear factor, defined as the ratio of the post- to preyield stiffness of the system; and z is the dimensionless Bouc–Wen hysteresis component.In testing processes, the actuator of the loading device was controlled by displacement input with a sinusoidal waveform as follows:(12)u=Apsin2πft,where Ap is the displacement amplitude; f is the excitation frequency; and t is the loading time.The main parameters of the MR damper are listed in Table5.Table 5
The main parameters of the MR damper.
Parameters
Values
Stroke (mm)
±
50
External diameter of cylinder (mm)
194
Internal diameter of cylinder (mm)
160
Diameter of the piston rod (mm)
80
Effective length of the piston (mm)
250
Coil (n×turns)
5
×
840
Height of the annular flow path (mm)
2
Trench diameter of coil (mm)
110
Coil resistance (Ω)
5
×
8.1In order to verify the effectiveness of the IMPSCO algorithm, four cases were simulated and discussed here. The details can be seen in Table6. Herein, the electricity value was set as1.2 A and the excitation frequency was 0.1 Hz.Table 6
Different cases for the MR damper experiment.
Case
Amplitude (mm) (Ap)
Case 1
5
Case 2
10
Case 3
20
Case 4
30


#### 3.3.2. Identification Process

As described in Section 3.1, the parameters to be identified should be encoded first. It is noted that the stiffness in the modified Bouc–Wen model (A) was redundant, demonstrated by Worden and Becker [28] using principled Bayesian approach, and thus A was set to 1. The search range of other corresponding parameters are listed in Table 7.Table 7
The parameters encoded in the IMPSCO for the MR damper experiment.
Parameters
Description
Initial range of parameters
f
0
Initial damping force
0,30000
c
Damping of MR damper
0,1.5
×
10
8
k
Stiffness of MR damper
0,3
×
10
8
α
Bilinear factor defined as the post- to preyield stiffness ratio of the system
0,3
×
10
8
β
The loop fatness of the system
0,30
γ
The loop pinching of the system
0,30
n
The abruptness parameter of the system
0,30Due to the fact that there was no acceleration record during the experiment, the error between measured displacement and produced displacement was set as the fitness function (equation (13)) and the parameters of IMPSCO in this study were the same as those in Section 3.2.2.(13)maxfθ=1∑i=1N∑kl=1LXitkl−Xitklp2,where Xt and Xtp are the reference measured displacement and the reproduced displacement, respectively.


#### 3.3.3. Identification Results and Discussion

According to the process described in Section3.2.2, we utilized the experimental displacement as described in Figure 7 to identify the parameters in the modified Bouc–Wen model. The identification results are presented in Table 8.Figure 7
Comparison between identified displacements and measured ones.Table 8
Identification results for the MR damper experiment.
Case
c (N·s/m)
k (N/m)
α
β
γ
N
f
0 (N)
Case 1
1.22 × 107
1.76 × 107
6.76 × 105
16.51
16.39
15.70
3286.92
Case 2
1.38 × 107
8.02 × 106
2.42 × 106
18.12
12.81
23.96
3668.17
Case 3
9.95 × 106
3.60 × 106
4.80 × 105
4.53
6.33
28.94
3000.00
Case 4
7.31 × 106
1.90 × 106
4.73 × 105
5.65
8.55
29.53
3000.00Due to the fact that it is difficult to measure the values of the parameters in the modified Bouc–Wen model matched with the MR damper, the identification results listed in Table8 were used to reconstruct the curve of displacements and damping forces so as to evaluate the efficiency of the proposed IMPSCO algorithm. The results are also shown in Figure 8. Herein, the identified values are also plotted in the same figure to compare with the measured ones.Figure 8
Comparison between identified values and measured values. (a)Ap=5mm. (b)Ap=10mm. (c)Ap=20mm. (d)Ap=30mm. (e) General damping force vs displacement curve.
(a)
(b)
(c)
(d)
(e)It can be seen from Figures7 and 8 that (1) the errors between identified values and measured ones for case 1 are much larger than those for cases 2, 3 and 4; (2) the identified values agree well with experimental values when the amplitude is larger than 5 mm, namely, cases 2, 3 and 4; and (3) the identified displacements coincide better with the measured values. The reason for the former is that when the amplitude of the MR damper is 5 mm, the damper just enter into the nonlinear phase. This implies that the modified Bouc–Wen model is not very suitable for simulating the MR damper in this early stage. By contrast, when the amplitude of the MR damper is larger than 5 mm and the MR damper comes into nonlinear phase, the identified results match well with experimental results. It is concluded that when the MR damper comes into the nonlinear phase, the modified Bouc–Wen model as described in equation (11) is perfect in depicting the force-displacement hysteretic relationship of the MR damper. In addition, the proposed IMPSCO algorithm is applicable to identifying parameters in the modified Bouc–Wen model.


## 3.1. Theory of the IMPSCO Algorithm

For basic multiparticle swarm coevolution optimization algorithm (MPSCO) [25], multiple subpopulations are divided into two layers. All particles from the upper layer follow the optimum of the entire population so as to obtain a faster convergence speed, while all particles from the lower layer follow the optimum of the subpopulation to ensure the population diversity. Although the performance of basic MPSCO is better than standard PSO in some respects, the subpopulations in the lower layer still perform the process of standard PSO, which makes particle falling into the local optimum possible. To solve this problem, we proposed the IMPSCO algorithm, and the details can be seen in Reference [21], and it was applied to locate and quantify damage in a linear structure. Figure 2 shows the flowchart of the IMPSCO algorithm.Figure 2
The flowchart of the IMPSCO algorithm.Wherein, the key point in IMPSCO is worst particle replacement. More specially, when the particle is recorded as the worst for the predetermined timesIw, replace it with local center of gravity Gg:(7)Gg=∑i=1szis,where s and zi represent the number of the selected excellent particles and their position, respectively.The most important task for optimization is to determine the fitness function in the IMPSCO. A useful fitness function not only has high accuracy but also needs less runtime, which are the main evaluation indices of the optimization algorithm for detecting sudden structural stiffness change. In addition, some parameters in the nonlinear hysteretic model may have little effects on the shape of the hysteretic curve. Therefore, we utilized the difference between the calculated and testing accelerations as the fitness function value, which contains lots of structural sensitive information. However, it is difficult to calculate the minimum. We employed the reciprocal of the mean square error (MSE) between the reproduced acceleration responsesX¨tp and the reference measured acceleration responses X¨t and replaced the mean square error (MSE) between them as the fitness function. When cast in discrete form, it can be expressed as(8)maxfθ=1∑i=1N∑kl=1LX¨itkl−X¨itklp2,where θ is the structural parameters vector to be optimized; N is the total number of measuring points; and L is the length of time history.


## 3.2. Feasibility: Numerical Simulation

### 3.2.1. Model Description

To evaluate the performance of IMPSCO algorithm applied to identification of nonlinear hysteretic base-isolation systems, a realistic base-isolated system is created as shown in Figure3 [26]. The numerical model can be simplified as a two degree-of-freedom (DOF) structure. Its basic definition includes the following: mass of each floor: m1=m2=1kg; initial stiffness of each floor: k1=k2=102.64N/m; damping coefficients of each floor: ζ1=ζ2=2%; parameters in the classical Bouc–Wen model: A=1;β=0.5;γ=0.5;n=2;α=0.1; the yield displacement of the system: Y=0.045m; and the basic period of structure: T=1s.Figure 3
2-DOF structure with Bouc–Wen model subjected to seismic excitation.Meanwhile, according to the sensitivity analysis of the parameters in the classical Bouc–Wen model listed in Reference [26], we assumed that some parameters with low sensitivity are constants in some cases in order to reduce the identified parameters in IMPSCO and simply the identification process. In summary, this study considered four cases with different parameters known. In addition, the search range of each parameter (θ) is within the interval of [0.5b, 2b], where b is the theoretical value of the corresponding parameter (θ).In this study, the structure was subjected to El-Centro earthquake wave, Kobe earthquake wave, and Northridge earthquake wave, respectively, at the base, as shown in Figure3. In order to reduce the particle dimension of IMPSCO, four cases were simulated to observe the feasibility of IMPSCO algorithm, as shown in Table 1. Therefore, by setting the sampling frequency of 50 Hz and the loading time t = 30 s, the structural acceleration responses were calculated by using the classical fourth-order Runge–Kutta algorithm [24].Table 1
Different cases for simulation.
Cases
Unknown parameters
Known parameters
S1
k
1
,
k
2
,
ζ
1
,
ζ
2
,
Y
,
α
,
β
,
γ
,
n
—
S2
k
1
,
k
2
,
ζ
1
,
ζ
2
,
Y
,
α
,
β
,
n
γ
=
0.5
S3
k
1
,
k
2
,
ζ
1
,
ζ
2
,
Y
,
β
,
n
α
=
0.1
;
γ
=
0.5
S4
k
1
,
k
2
,
ζ
1
,
ζ
2
,
Y
,
α
,
β
,
γ
n
=
2To consider the operating environment, Gaussian white noise with different levels was added to all structural responses and the earthquake motion. The structural responses were then represented as(9)yi=yia×1+εR,where yi and yia represent the contaminated and theoretical signals, respectively; R is a normally distributed random variable with zero mean and a derivation of 1; and ε is an index representing the noise level.Therefore, the signal-to-noise ratio (SNR) is defined as(10)SNR=20log101ε=20log10AsignalAnoise,where Asignal and Anoise are the amplitudes of the signal and noise, respectively.In this paper, the SNR was 30 dB.


### 3.2.2. Identification Process

In order to identify parameters of 2-DOF structure, two issues need to be ensured for the IMPSCO algorithm: encoding parameters and initial setting parameters. The details can be seen as follows.(1) Parameter Encoding Involved in the IMPSCO. This study assumed that the mass of each floor was known. Therefore, the parameters encoded in the IMPSCO were initial stiffness (k1, k2) and damping coefficients (ζ1, ζ2), yield displacement (Y) of the first floor, bilinear factor (α), loop fatness (β), loop pinching (γ), and abruptness parameters (n) in the classical Bouc–Wen model system, and the details could be seen in Figure 4. Consequently, the number of subpopulations m was 9.Figure 4
The parameters encoded in the IMPSCO.(2) Parameters in IMPSCO. To maximize the fitness function using IMPSCO presented in Section 3.1, some parameters should be set and initialized in advance for IMPSCO. Generally, 3∼5 subpopulations with 30 to 100 particles in total were sufficient. The inertia weight w was set linearly varying from 0.9 to 0.4 and learning factors c1 and c2 were set to 2 simultaneously. In addition, the limited times for the worst record Iw and the number of the selected excellent particles s were both empirically set to 5∼10, respectively. Finally, the parameters of IMPSCO in this study were initialized as shown in Table 2.Table 2
Specifications of IMPSCO algorithm.
Parameter
Description
Value
mm
The total number of subpopulations
3
nn
The size of a subpopulation
10
c
1
Cognitive parameter
2
c
2
Social parameter
2
ρ
min
Minimum of inertia weight
0.4
ρ
max
Maximum of inertia weight
0.9
I
w
Limited times for worst record
5
S
The number of the selected excellent particles
6
I
max
The maximum iteration times
200


### 3.2.3. Identification Results and Discussion

According to the process described in Section3.1, the responses shown in Figure 5 were employed to identify the parameters in the nonlinear structure and in Bouc–Wen model system. The results are listed in Table 3, and relative errors between identified results and theoretical values are shown in Table 4.Figure 5
Acceleration responses of the simulated structure under different earthquake waves,∗a1 is the acceleration responses of the first floor; a2 is the acceleration responses of the second floor; aq is corresponding earthquake excitation. (a) El-Centro earthquake waves. (b) Kobe earthquake waves. (c) Northridge earthquake waves.
(a)
(b)
(c)Table 3
Identification results of the simulated structure subjected to the different earthquake excitation in different cases.
Case
Seismic excitation
k
1 (N/m)
k
2 (N/m)
ζ
1 (%)
ζ
2 (%)
Y (mm)
α
β
γ
n
S1
(a) El-Centro earthquake
102.45
102.73
1.86
2.06
46.46
0.09
0.55
0.52
2.06
(b) Kobe earthquake
102.79
102.73
2.12
1.83
47.38
0.10
0.56
0.57
1.95
(c) Northridge earthquake
103.01
102.58
1.96
1.98
49.60
0.10
0.61
0.61
1.97
S2
(a) El-Centro earthquake
103.15
102.46
2.09
2.16
44.13
0.11
0.47
—
2.01
(b) Kobe earthquake
102.65
102.74
2.12
1.82
44.65
0.10
0.49
—
1.96
(c) Northridge earthquake
101.46
103.09
2.38
2.09
45.19
0.10
0.46
—
2.07
S3
(a) El-Centro earthquake
102.86
102.70
2.03
2.24
44.70
—
0.48
—
2.01
(b) Kobe earthquake
103.13
102.71
2.15
1.84
44.39
—
0.48
—
1.93
(c) Northridge earthquake
102.72
103.01
2.28
2.05
44.71
—
0.46
—
1.97
S4
(a) El-Centro earthquake
103.06
102.55
2.08
2.17
47.78
0.11
0.55
0.61
—
(b) Kobe earthquake
102.54
102.73
2.16
1.84
47.62
0.10
0.56
0.58
—
(c) Northridge earthquake
102.94
102.65
2.16
2.01
46.23
0.10
0.52
0.56
—Table 4
The relative errors between identified parameters and actual model parameters.
Case
Seismic excitation
Relative error (%)
k
1
k
2
ζ
1
ζ
2
Y
α
β
γ
n
S1
(a) El-Centro earthquake
0.19
0.09
7.00
3.00
3.24
10.00
10.00
4.00
3.00
(b) Kobe earthquake
0.15
0.09
6.00
8.50
5.29
0.00
12.00
14.00
2.50
(c) Northridge earthquake
0.36
0.06
2.00
1.00
10.22
0.00
22.00
22.00
1.50
Average values
0.23
0.08
5.00
4.17
6.25
3.33
14.67
13.33
2.33
S2
(a) El-Centro earthquake
0.50
0.18
4.50
8.00
1.93
10.00
6.00
—
0.50
(b) Kobe earthquake
0.01
0.10
6.00
9.00
0.78
0.00
2.00
—
2.00
(c) Northridge earthquake
1.15
0.44
19.00
4.50
0.42
0.00
8.00
—
3.50
Average values
0.55
0.24
9.83
7.17
1.04
3.33
5.33
—
2.00
S3
(a) El-Centro earthquake
0.21
0.06
1.50
12.00
0.67
—
4.00
—
0.50
(b) Kobe earthquake
0.48
0.07
7.50
8.00
1.36
—
4.00
—
3.50
(c) Northridge earthquake
0.08
0.36
14.00
2.50
0.64
—
8.00
—
1.50
Average values
0.26
0.16
7.67
7.50
0.89
—
5.33
—
1.83
S4
(a) El-Centro earthquake
0.41
0.09
4.00
8.50
6.18
10.00
10.00
22.00
—
(b) Kobe earthquake
0.10
0.09
8.00
8.00
5.82
0.00
12.00
16.00
—
(c) Northridge earthquake
0.29
0.01
8.00
0.50
2.73
0.00
4.00
12.00
—
Average values
0.27
0.06
6.67
5.67
4.91
3.33
8.67
16.67
—The following conclusions and remarks can be made from the two tables:(1)
The relative errors of initial stiffness of system (k) and the abruptness parameter (n) in the classical Bouc–Wen model are less than 3.5%, and the rest of identified parameters have lower computational accuracy. More specifically, the maximum relative errors for stiffness (α), loop fatness (β), and loop pinching (γ) are 10.00% (S1, S2, S4), 22.00% (S1), and 22.00% (S1, S4), respectively; the corresponding values are 10.22% (S1), 1.93% (S2), 1.36% (S3), and 6.18% (S4) for the yield displacement (Y); in addition, the maximum relative errors for damping coefficients (ζ1 and ζ2) are 8.50% (S1), 19.00% (S2), 14.00% (S3), and 8.50% (S4), respectively.There are two reasons for this phenomenon. The first one is that when we used the error between reproduced and measured acceleration responses as the fitness function, the stiffness of the system mainly influences the reproduced accelerations. By contrast, the damping coefficients have little effect on the reproduced accelerations. As a result, it is difficult to accurately identify the damping coefficients. The second reason is that, as presented in Reference [27], the abruptness parameter in the classical Bouc–Wen model (n) is a power factor, which changes the smoothness of transition from elastic region to plastic region in the hysteretic versus actual displacement curves, with little change in the acceleration response. Correspondingly, the remaining three parameters in the classical Bouc–Wen model are not sensitive to accelerations.(2)
It is accurate to identify stiffness of the system (k1,k2) and the abruptness parameter (n) in the classical Bouc–Wen model by using IMPSCO with the parameter listed in Table 2. Furthermore, for four cases (with different parameters to be identified and different particle dimension of IMPSCO), all the parameters of the nonlinear structure can be basically identified by IMPSCO. Although the identified results for the parameters in the classical Bouc–Wen model have lower accuracy than that of stiffness, the IMPSCO algorithm is suitable for identification of sudden stiffness change in actual structure. Just as shown in Table 4, the overall ranking for stiffness (k1,k2) and damping coefficients (ζ1,ζ2) is S1 > S4 > S3 > S2, and the overall ranking for remaining parameters is S3 > S2 > S4 > S1. This indicates that for all the parameters unknown case (S1), the identified results for the stiffness and damping have excellent accuracy. By contrast, the identified results for the parameters in the Bouc–Wen model have lower accuracy than other cases. However, the maximum error in the identification process is less than 16.7%. These error values are well within construction errors. It also demonstrates that the IMPSCO algorithm is applicable for detecting structural damage induced by stiffness change in structures.


## 3.2.1. Model Description

To evaluate the performance of IMPSCO algorithm applied to identification of nonlinear hysteretic base-isolation systems, a realistic base-isolated system is created as shown in Figure3 [26]. The numerical model can be simplified as a two degree-of-freedom (DOF) structure. Its basic definition includes the following: mass of each floor: m1=m2=1kg; initial stiffness of each floor: k1=k2=102.64N/m; damping coefficients of each floor: ζ1=ζ2=2%; parameters in the classical Bouc–Wen model: A=1;β=0.5;γ=0.5;n=2;α=0.1; the yield displacement of the system: Y=0.045m; and the basic period of structure: T=1s.Figure 3
2-DOF structure with Bouc–Wen model subjected to seismic excitation.Meanwhile, according to the sensitivity analysis of the parameters in the classical Bouc–Wen model listed in Reference [26], we assumed that some parameters with low sensitivity are constants in some cases in order to reduce the identified parameters in IMPSCO and simply the identification process. In summary, this study considered four cases with different parameters known. In addition, the search range of each parameter (θ) is within the interval of [0.5b, 2b], where b is the theoretical value of the corresponding parameter (θ).In this study, the structure was subjected to El-Centro earthquake wave, Kobe earthquake wave, and Northridge earthquake wave, respectively, at the base, as shown in Figure3. In order to reduce the particle dimension of IMPSCO, four cases were simulated to observe the feasibility of IMPSCO algorithm, as shown in Table 1. Therefore, by setting the sampling frequency of 50 Hz and the loading time t = 30 s, the structural acceleration responses were calculated by using the classical fourth-order Runge–Kutta algorithm [24].Table 1
Different cases for simulation.
Cases
Unknown parameters
Known parameters
S1
k
1
,
k
2
,
ζ
1
,
ζ
2
,
Y
,
α
,
β
,
γ
,
n
—
S2
k
1
,
k
2
,
ζ
1
,
ζ
2
,
Y
,
α
,
β
,
n
γ
=
0.5
S3
k
1
,
k
2
,
ζ
1
,
ζ
2
,
Y
,
β
,
n
α
=
0.1
;
γ
=
0.5
S4
k
1
,
k
2
,
ζ
1
,
ζ
2
,
Y
,
α
,
β
,
γ
n
=
2To consider the operating environment, Gaussian white noise with different levels was added to all structural responses and the earthquake motion. The structural responses were then represented as(9)yi=yia×1+εR,where yi and yia represent the contaminated and theoretical signals, respectively; R is a normally distributed random variable with zero mean and a derivation of 1; and ε is an index representing the noise level.Therefore, the signal-to-noise ratio (SNR) is defined as(10)SNR=20log101ε=20log10AsignalAnoise,where Asignal and Anoise are the amplitudes of the signal and noise, respectively.In this paper, the SNR was 30 dB.


## 3.2.2. Identification Process

In order to identify parameters of 2-DOF structure, two issues need to be ensured for the IMPSCO algorithm: encoding parameters and initial setting parameters. The details can be seen as follows.(1) Parameter Encoding Involved in the IMPSCO. This study assumed that the mass of each floor was known. Therefore, the parameters encoded in the IMPSCO were initial stiffness (k1, k2) and damping coefficients (ζ1, ζ2), yield displacement (Y) of the first floor, bilinear factor (α), loop fatness (β), loop pinching (γ), and abruptness parameters (n) in the classical Bouc–Wen model system, and the details could be seen in Figure 4. Consequently, the number of subpopulations m was 9.Figure 4
The parameters encoded in the IMPSCO.(2) Parameters in IMPSCO. To maximize the fitness function using IMPSCO presented in Section 3.1, some parameters should be set and initialized in advance for IMPSCO. Generally, 3∼5 subpopulations with 30 to 100 particles in total were sufficient. The inertia weight w was set linearly varying from 0.9 to 0.4 and learning factors c1 and c2 were set to 2 simultaneously. In addition, the limited times for the worst record Iw and the number of the selected excellent particles s were both empirically set to 5∼10, respectively. Finally, the parameters of IMPSCO in this study were initialized as shown in Table 2.Table 2
Specifications of IMPSCO algorithm.
Parameter
Description
Value
mm
The total number of subpopulations
3
nn
The size of a subpopulation
10
c
1
Cognitive parameter
2
c
2
Social parameter
2
ρ
min
Minimum of inertia weight
0.4
ρ
max
Maximum of inertia weight
0.9
I
w
Limited times for worst record
5
S
The number of the selected excellent particles
6
I
max
The maximum iteration times
200


## 3.2.3. Identification Results and Discussion

According to the process described in Section3.1, the responses shown in Figure 5 were employed to identify the parameters in the nonlinear structure and in Bouc–Wen model system. The results are listed in Table 3, and relative errors between identified results and theoretical values are shown in Table 4.Figure 5
Acceleration responses of the simulated structure under different earthquake waves,∗a1 is the acceleration responses of the first floor; a2 is the acceleration responses of the second floor; aq is corresponding earthquake excitation. (a) El-Centro earthquake waves. (b) Kobe earthquake waves. (c) Northridge earthquake waves.
(a)
(b)
(c)Table 3
Identification results of the simulated structure subjected to the different earthquake excitation in different cases.
Case
Seismic excitation
k
1 (N/m)
k
2 (N/m)
ζ
1 (%)
ζ
2 (%)
Y (mm)
α
β
γ
n
S1
(a) El-Centro earthquake
102.45
102.73
1.86
2.06
46.46
0.09
0.55
0.52
2.06
(b) Kobe earthquake
102.79
102.73
2.12
1.83
47.38
0.10
0.56
0.57
1.95
(c) Northridge earthquake
103.01
102.58
1.96
1.98
49.60
0.10
0.61
0.61
1.97
S2
(a) El-Centro earthquake
103.15
102.46
2.09
2.16
44.13
0.11
0.47
—
2.01
(b) Kobe earthquake
102.65
102.74
2.12
1.82
44.65
0.10
0.49
—
1.96
(c) Northridge earthquake
101.46
103.09
2.38
2.09
45.19
0.10
0.46
—
2.07
S3
(a) El-Centro earthquake
102.86
102.70
2.03
2.24
44.70
—
0.48
—
2.01
(b) Kobe earthquake
103.13
102.71
2.15
1.84
44.39
—
0.48
—
1.93
(c) Northridge earthquake
102.72
103.01
2.28
2.05
44.71
—
0.46
—
1.97
S4
(a) El-Centro earthquake
103.06
102.55
2.08
2.17
47.78
0.11
0.55
0.61
—
(b) Kobe earthquake
102.54
102.73
2.16
1.84
47.62
0.10
0.56
0.58
—
(c) Northridge earthquake
102.94
102.65
2.16
2.01
46.23
0.10
0.52
0.56
—Table 4
The relative errors between identified parameters and actual model parameters.
Case
Seismic excitation
Relative error (%)
k
1
k
2
ζ
1
ζ
2
Y
α
β
γ
n
S1
(a) El-Centro earthquake
0.19
0.09
7.00
3.00
3.24
10.00
10.00
4.00
3.00
(b) Kobe earthquake
0.15
0.09
6.00
8.50
5.29
0.00
12.00
14.00
2.50
(c) Northridge earthquake
0.36
0.06
2.00
1.00
10.22
0.00
22.00
22.00
1.50
Average values
0.23
0.08
5.00
4.17
6.25
3.33
14.67
13.33
2.33
S2
(a) El-Centro earthquake
0.50
0.18
4.50
8.00
1.93
10.00
6.00
—
0.50
(b) Kobe earthquake
0.01
0.10
6.00
9.00
0.78
0.00
2.00
—
2.00
(c) Northridge earthquake
1.15
0.44
19.00
4.50
0.42
0.00
8.00
—
3.50
Average values
0.55
0.24
9.83
7.17
1.04
3.33
5.33
—
2.00
S3
(a) El-Centro earthquake
0.21
0.06
1.50
12.00
0.67
—
4.00
—
0.50
(b) Kobe earthquake
0.48
0.07
7.50
8.00
1.36
—
4.00
—
3.50
(c) Northridge earthquake
0.08
0.36
14.00
2.50
0.64
—
8.00
—
1.50
Average values
0.26
0.16
7.67
7.50
0.89
—
5.33
—
1.83
S4
(a) El-Centro earthquake
0.41
0.09
4.00
8.50
6.18
10.00
10.00
22.00
—
(b) Kobe earthquake
0.10
0.09
8.00
8.00
5.82
0.00
12.00
16.00
—
(c) Northridge earthquake
0.29
0.01
8.00
0.50
2.73
0.00
4.00
12.00
—
Average values
0.27
0.06
6.67
5.67
4.91
3.33
8.67
16.67
—The following conclusions and remarks can be made from the two tables:(1)
The relative errors of initial stiffness of system (k) and the abruptness parameter (n) in the classical Bouc–Wen model are less than 3.5%, and the rest of identified parameters have lower computational accuracy. More specifically, the maximum relative errors for stiffness (α), loop fatness (β), and loop pinching (γ) are 10.00% (S1, S2, S4), 22.00% (S1), and 22.00% (S1, S4), respectively; the corresponding values are 10.22% (S1), 1.93% (S2), 1.36% (S3), and 6.18% (S4) for the yield displacement (Y); in addition, the maximum relative errors for damping coefficients (ζ1 and ζ2) are 8.50% (S1), 19.00% (S2), 14.00% (S3), and 8.50% (S4), respectively.There are two reasons for this phenomenon. The first one is that when we used the error between reproduced and measured acceleration responses as the fitness function, the stiffness of the system mainly influences the reproduced accelerations. By contrast, the damping coefficients have little effect on the reproduced accelerations. As a result, it is difficult to accurately identify the damping coefficients. The second reason is that, as presented in Reference [27], the abruptness parameter in the classical Bouc–Wen model (n) is a power factor, which changes the smoothness of transition from elastic region to plastic region in the hysteretic versus actual displacement curves, with little change in the acceleration response. Correspondingly, the remaining three parameters in the classical Bouc–Wen model are not sensitive to accelerations.(2)
It is accurate to identify stiffness of the system (k1,k2) and the abruptness parameter (n) in the classical Bouc–Wen model by using IMPSCO with the parameter listed in Table 2. Furthermore, for four cases (with different parameters to be identified and different particle dimension of IMPSCO), all the parameters of the nonlinear structure can be basically identified by IMPSCO. Although the identified results for the parameters in the classical Bouc–Wen model have lower accuracy than that of stiffness, the IMPSCO algorithm is suitable for identification of sudden stiffness change in actual structure. Just as shown in Table 4, the overall ranking for stiffness (k1,k2) and damping coefficients (ζ1,ζ2) is S1 > S4 > S3 > S2, and the overall ranking for remaining parameters is S3 > S2 > S4 > S1. This indicates that for all the parameters unknown case (S1), the identified results for the stiffness and damping have excellent accuracy. By contrast, the identified results for the parameters in the Bouc–Wen model have lower accuracy than other cases. However, the maximum error in the identification process is less than 16.7%. These error values are well within construction errors. It also demonstrates that the IMPSCO algorithm is applicable for detecting structural damage induced by stiffness change in structures.


## 3.3. Feasibility: Experiment

### 3.3.1. Test Setup

A MR damper test [27] was presented to further validate the proposed IMPSCO algorithm. The test’s configuration schematic and the photograph of the testing system are shown in Figure 6.Figure 6
MR damper performance test setup [27].It is noted that the modified Bouc–Wen model is superior in describing the force-displacement relationship of MR damper. Its expression is as follows:(11)F=cx˙+kx+αz+f0,z˙=Ax˙−βx˙zzn−1+γx˙zn,where F and f0 are the damping force and initial damping force; c and k are damping and the stiffness of the system, respectively; α is the bilinear factor, defined as the ratio of the post- to preyield stiffness of the system; and z is the dimensionless Bouc–Wen hysteresis component.In testing processes, the actuator of the loading device was controlled by displacement input with a sinusoidal waveform as follows:(12)u=Apsin2πft,where Ap is the displacement amplitude; f is the excitation frequency; and t is the loading time.The main parameters of the MR damper are listed in Table5.Table 5
The main parameters of the MR damper.
Parameters
Values
Stroke (mm)
±
50
External diameter of cylinder (mm)
194
Internal diameter of cylinder (mm)
160
Diameter of the piston rod (mm)
80
Effective length of the piston (mm)
250
Coil (n×turns)
5
×
840
Height of the annular flow path (mm)
2
Trench diameter of coil (mm)
110
Coil resistance (Ω)
5
×
8.1In order to verify the effectiveness of the IMPSCO algorithm, four cases were simulated and discussed here. The details can be seen in Table6. Herein, the electricity value was set as1.2 A and the excitation frequency was 0.1 Hz.Table 6
Different cases for the MR damper experiment.
Case
Amplitude (mm) (Ap)
Case 1
5
Case 2
10
Case 3
20
Case 4
30


### 3.3.2. Identification Process

As described in Section 3.1, the parameters to be identified should be encoded first. It is noted that the stiffness in the modified Bouc–Wen model (A) was redundant, demonstrated by Worden and Becker [28] using principled Bayesian approach, and thus A was set to 1. The search range of other corresponding parameters are listed in Table 7.Table 7
The parameters encoded in the IMPSCO for the MR damper experiment.
Parameters
Description
Initial range of parameters
f
0
Initial damping force
0,30000
c
Damping of MR damper
0,1.5
×
10
8
k
Stiffness of MR damper
0,3
×
10
8
α
Bilinear factor defined as the post- to preyield stiffness ratio of the system
0,3
×
10
8
β
The loop fatness of the system
0,30
γ
The loop pinching of the system
0,30
n
The abruptness parameter of the system
0,30Due to the fact that there was no acceleration record during the experiment, the error between measured displacement and produced displacement was set as the fitness function (equation (13)) and the parameters of IMPSCO in this study were the same as those in Section 3.2.2.(13)maxfθ=1∑i=1N∑kl=1LXitkl−Xitklp2,where Xt and Xtp are the reference measured displacement and the reproduced displacement, respectively.


### 3.3.3. Identification Results and Discussion

According to the process described in Section3.2.2, we utilized the experimental displacement as described in Figure 7 to identify the parameters in the modified Bouc–Wen model. The identification results are presented in Table 8.Figure 7
Comparison between identified displacements and measured ones.Table 8
Identification results for the MR damper experiment.
Case
c (N·s/m)
k (N/m)
α
β
γ
N
f
0 (N)
Case 1
1.22 × 107
1.76 × 107
6.76 × 105
16.51
16.39
15.70
3286.92
Case 2
1.38 × 107
8.02 × 106
2.42 × 106
18.12
12.81
23.96
3668.17
Case 3
9.95 × 106
3.60 × 106
4.80 × 105
4.53
6.33
28.94
3000.00
Case 4
7.31 × 106
1.90 × 106
4.73 × 105
5.65
8.55
29.53
3000.00Due to the fact that it is difficult to measure the values of the parameters in the modified Bouc–Wen model matched with the MR damper, the identification results listed in Table8 were used to reconstruct the curve of displacements and damping forces so as to evaluate the efficiency of the proposed IMPSCO algorithm. The results are also shown in Figure 8. Herein, the identified values are also plotted in the same figure to compare with the measured ones.Figure 8
Comparison between identified values and measured values. (a)Ap=5mm. (b)Ap=10mm. (c)Ap=20mm. (d)Ap=30mm. (e) General damping force vs displacement curve.
(a)
(b)
(c)
(d)
(e)It can be seen from Figures7 and 8 that (1) the errors between identified values and measured ones for case 1 are much larger than those for cases 2, 3 and 4; (2) the identified values agree well with experimental values when the amplitude is larger than 5 mm, namely, cases 2, 3 and 4; and (3) the identified displacements coincide better with the measured values. The reason for the former is that when the amplitude of the MR damper is 5 mm, the damper just enter into the nonlinear phase. This implies that the modified Bouc–Wen model is not very suitable for simulating the MR damper in this early stage. By contrast, when the amplitude of the MR damper is larger than 5 mm and the MR damper comes into nonlinear phase, the identified results match well with experimental results. It is concluded that when the MR damper comes into the nonlinear phase, the modified Bouc–Wen model as described in equation (11) is perfect in depicting the force-displacement hysteretic relationship of the MR damper. In addition, the proposed IMPSCO algorithm is applicable to identifying parameters in the modified Bouc–Wen model.


## 3.3.1. Test Setup

A MR damper test [27] was presented to further validate the proposed IMPSCO algorithm. The test’s configuration schematic and the photograph of the testing system are shown in Figure 6.Figure 6
MR damper performance test setup [27].It is noted that the modified Bouc–Wen model is superior in describing the force-displacement relationship of MR damper. Its expression is as follows:(11)F=cx˙+kx+αz+f0,z˙=Ax˙−βx˙zzn−1+γx˙zn,where F and f0 are the damping force and initial damping force; c and k are damping and the stiffness of the system, respectively; α is the bilinear factor, defined as the ratio of the post- to preyield stiffness of the system; and z is the dimensionless Bouc–Wen hysteresis component.In testing processes, the actuator of the loading device was controlled by displacement input with a sinusoidal waveform as follows:(12)u=Apsin2πft,where Ap is the displacement amplitude; f is the excitation frequency; and t is the loading time.The main parameters of the MR damper are listed in Table5.Table 5
The main parameters of the MR damper.
Parameters
Values
Stroke (mm)
±
50
External diameter of cylinder (mm)
194
Internal diameter of cylinder (mm)
160
Diameter of the piston rod (mm)
80
Effective length of the piston (mm)
250
Coil (n×turns)
5
×
840
Height of the annular flow path (mm)
2
Trench diameter of coil (mm)
110
Coil resistance (Ω)
5
×
8.1In order to verify the effectiveness of the IMPSCO algorithm, four cases were simulated and discussed here. The details can be seen in Table6. Herein, the electricity value was set as1.2 A and the excitation frequency was 0.1 Hz.Table 6
Different cases for the MR damper experiment.
Case
Amplitude (mm) (Ap)
Case 1
5
Case 2
10
Case 3
20
Case 4
30


## 3.3.2. Identification Process

As described in Section 3.1, the parameters to be identified should be encoded first. It is noted that the stiffness in the modified Bouc–Wen model (A) was redundant, demonstrated by Worden and Becker [28] using principled Bayesian approach, and thus A was set to 1. The search range of other corresponding parameters are listed in Table 7.Table 7
The parameters encoded in the IMPSCO for the MR damper experiment.
Parameters
Description
Initial range of parameters
f
0
Initial damping force
0,30000
c
Damping of MR damper
0,1.5
×
10
8
k
Stiffness of MR damper
0,3
×
10
8
α
Bilinear factor defined as the post- to preyield stiffness ratio of the system
0,3
×
10
8
β
The loop fatness of the system
0,30
γ
The loop pinching of the system
0,30
n
The abruptness parameter of the system
0,30Due to the fact that there was no acceleration record during the experiment, the error between measured displacement and produced displacement was set as the fitness function (equation (13)) and the parameters of IMPSCO in this study were the same as those in Section 3.2.2.(13)maxfθ=1∑i=1N∑kl=1LXitkl−Xitklp2,where Xt and Xtp are the reference measured displacement and the reproduced displacement, respectively.


## 3.3.3. Identification Results and Discussion

According to the process described in Section3.2.2, we utilized the experimental displacement as described in Figure 7 to identify the parameters in the modified Bouc–Wen model. The identification results are presented in Table 8.Figure 7
Comparison between identified displacements and measured ones.Table 8
Identification results for the MR damper experiment.
Case
c (N·s/m)
k (N/m)
α
β
γ
N
f
0 (N)
Case 1
1.22 × 107
1.76 × 107
6.76 × 105
16.51
16.39
15.70
3286.92
Case 2
1.38 × 107
8.02 × 106
2.42 × 106
18.12
12.81
23.96
3668.17
Case 3
9.95 × 106
3.60 × 106
4.80 × 105
4.53
6.33
28.94
3000.00
Case 4
7.31 × 106
1.90 × 106
4.73 × 105
5.65
8.55
29.53
3000.00Due to the fact that it is difficult to measure the values of the parameters in the modified Bouc–Wen model matched with the MR damper, the identification results listed in Table8 were used to reconstruct the curve of displacements and damping forces so as to evaluate the efficiency of the proposed IMPSCO algorithm. The results are also shown in Figure 8. Herein, the identified values are also plotted in the same figure to compare with the measured ones.Figure 8
Comparison between identified values and measured values. (a)Ap=5mm. (b)Ap=10mm. (c)Ap=20mm. (d)Ap=30mm. (e) General damping force vs displacement curve.
(a)
(b)
(c)
(d)
(e)It can be seen from Figures7 and 8 that (1) the errors between identified values and measured ones for case 1 are much larger than those for cases 2, 3 and 4; (2) the identified values agree well with experimental values when the amplitude is larger than 5 mm, namely, cases 2, 3 and 4; and (3) the identified displacements coincide better with the measured values. The reason for the former is that when the amplitude of the MR damper is 5 mm, the damper just enter into the nonlinear phase. This implies that the modified Bouc–Wen model is not very suitable for simulating the MR damper in this early stage. By contrast, when the amplitude of the MR damper is larger than 5 mm and the MR damper comes into nonlinear phase, the identified results match well with experimental results. It is concluded that when the MR damper comes into the nonlinear phase, the modified Bouc–Wen model as described in equation (11) is perfect in depicting the force-displacement hysteretic relationship of the MR damper. In addition, the proposed IMPSCO algorithm is applicable to identifying parameters in the modified Bouc–Wen model.


## 4. Two-Stage Identification Strategy for Sudden Stiffness Change in Nonlinear Hysteretic Structures

The improved MPSCO algorithm is able to identify structural stiffness, but it cannot detect when the structure stiffness changed. As a result, a new two-stage damage detection strategy by integrating DWT-FastICA and the IMPSCO algorithm is proposed to identify the sudden stiffness change of the acceleration response in a nonlinear hysteric structure. The schematic diagram of the damage detection strategy is depicted in Figure9. In the first stage, the measured structural dynamic responses are preprocessed firstly by discrete wavelet transform (DWT), and then the fast independent component analysis (Fast ICA) is employed to detect the structural response abnormality and initially locate damage location. In the second stage, the structural responses are divided into segments which are used to identify the time-varying nonlinear structural parameters by the IMPSCO algorithm, whereby the damage location and extent can be identified accurately.Figure 9
Two-stage sudden stiffness change of nonlinear structure detection strategy. (a) The first stage: primary damage assessment. (b) The second stage: refined damage assessment.
(a)
(b)


### 4.1. The First Stage: Primary Damage Assessment

Two essential tasks in the first stage are detecting damage novelty and initially locating damage.


#### 4.1.1. Separation of Wavelet Detailed Components

To complete the damage novelty detection, the time-series dynamic responses is preprocessed by DWT in order to separate the wavelet detailed components. DWT [29] is a useful tool for time-frequency analysis, which achieves multiresolution analysis of a signal by decomposing it into high-frequency (detail) and low-frequency (approximation) components at each level. If a signal is decomposed into LL levels, it can be reconstructed by(14)X¨t=∑ll=1LLX¨dllt+X¨aLLt=∑ll=1LL∑kk=−∞∞dkk,llψkk,llt+∑kk=−∞∞ckk,LLϕkk,LLt,where X¨dllt is the detailed component at the ll-th scale level; X¨aLLt is the approximated component at the LL-th scale level; kk and ll are the translation and scale parameters, respectively; dkk,ll and ψkk,llt are the detail coefficient and wavelet function at the ll-th scale level, respectively; and ckl,LL and ϕkk,LLt are the approximation coefficient and scaling function at the LL-th scale level, respectively.Two parameters should be predetermined using DWT to detect the singularity (or discontinuity points) of a signal. The first one is selecting the proper wavelet basis and the other one is the decomposition level involving in DWT. After determining the parameters, DWT is, respectively, performed in the structural dynamic responsesXi (i = 1, 2, …, m; m is the total number of sensors) measured from each sensor location. Consequently, the wavelet detailed components wXi on a certain scale which contain the most damage information can be separated, respectively. So, the wavelet detailed signal matrix can be built as wX = [wX1, wX2 , …, wXm]T.


#### 4.1.2. Extraction of Feature Independent Component

A new singularity analysis method called DWT-ICA by integrating DWT with independent component analysis (ICA) was recently proposed and applied to structural damage detection. This method shows excellent performance in identifying sudden stiffness change for linear structures [6]. Independent component analysis (ICA) [30] is a computational technique for feature extraction in signal processing, which can separate a multivariate observed signal to statistically independent sources. It can be specially described as(15)wXt=AwSt,where wSt=wS1t,…,wSmt represents an n-dimensional independent source matrix and A represents an m × n (m ≥ n) mixing matrix.The principle of ICA can be summarized as an optimization process to search for proper estimation of the inverse ofA (i.e., demixing matrix W) such that the source signal matrix wS can be recovered by(16)wSt=WwXt.Therefore, there are two key issues in the ICA: (1) the proper objective function in order to decide whether the obtained source signals are statistically independent; (2) the effective algorithm to implement the optimization of the objective function.The Fast ICA proposed by Hyvärinen and Oja [31] is one of the most efficient and popular algorithms in revised ICA. It is based on a fixed-point iteration scheme maximizing non-Gaussianity as the objective function.To extract the damage information hidden in the noisy wavelet-domain signals, the detailed signal matrixwX is processed by using the Fast ICA, thereby nIC independent components ICj (j = 1, 2, …, nIC) can be obtained. Subsequently, the feature independent component (FIC) which contains the spike property can be extracted, and the position of the spike in the time history indicates the damage occurrence instant.


#### 4.1.3. Initial Localization of Damage

After conducting the Fast ICA, the wavelet detailed signal matrixwX can be reconstructed by(17)wX=A×IC=∑i=1naj⋅ICj,where aj = [a1j, …, aij, …, amj]T is the mixing vector which is called as the source distribution vector (SDV). Herein, the element aij that is defined as the source distribution factor (SDF) represents the proportion of ICj distributed in wXi; ICj is the j-th independent component extracted by the Fast ICA.When ICj contained structural singularity information, it is called the feature independent component (FICk k = 1, 2, …, h). If aij is the largest among its SDV, it indicates its corresponding location, where the response measured in the i-th sensor location Xi contains the most FIC, is the possibly damaged location.The SDV can be normalized further by equation (18). Thus, the structural damage can be initially localized by the maximum among the NSDV of the FIC.(18)NSDF=NSDF1,NSDF2,⋯,NSDFi,NSDFmT,NSDFi=aij∑i=1maij,i=1,2,…,m.


### 4.2. Refined Damage Assessment

First, the original time-series responses are divided intoks + 1 segments as X¨t1,X¨t2,…,X¨tks+1 according to the damage instant identified through FICks, ks = 1, 2, ..., h.Second, IMPSCO is performed on theks + 1 segment X¨t1,X¨t2,…,X¨tks+1, respectively; thus, the damage location and extent can be identified precisely according to the stiffness variation at each time period.


## 4.1. The First Stage: Primary Damage Assessment

Two essential tasks in the first stage are detecting damage novelty and initially locating damage.


### 4.1.1. Separation of Wavelet Detailed Components

To complete the damage novelty detection, the time-series dynamic responses is preprocessed by DWT in order to separate the wavelet detailed components. DWT [29] is a useful tool for time-frequency analysis, which achieves multiresolution analysis of a signal by decomposing it into high-frequency (detail) and low-frequency (approximation) components at each level. If a signal is decomposed into LL levels, it can be reconstructed by(14)X¨t=∑ll=1LLX¨dllt+X¨aLLt=∑ll=1LL∑kk=−∞∞dkk,llψkk,llt+∑kk=−∞∞ckk,LLϕkk,LLt,where X¨dllt is the detailed component at the ll-th scale level; X¨aLLt is the approximated component at the LL-th scale level; kk and ll are the translation and scale parameters, respectively; dkk,ll and ψkk,llt are the detail coefficient and wavelet function at the ll-th scale level, respectively; and ckl,LL and ϕkk,LLt are the approximation coefficient and scaling function at the LL-th scale level, respectively.Two parameters should be predetermined using DWT to detect the singularity (or discontinuity points) of a signal. The first one is selecting the proper wavelet basis and the other one is the decomposition level involving in DWT. After determining the parameters, DWT is, respectively, performed in the structural dynamic responsesXi (i = 1, 2, …, m; m is the total number of sensors) measured from each sensor location. Consequently, the wavelet detailed components wXi on a certain scale which contain the most damage information can be separated, respectively. So, the wavelet detailed signal matrix can be built as wX = [wX1, wX2 , …, wXm]T.


### 4.1.2. Extraction of Feature Independent Component

A new singularity analysis method called DWT-ICA by integrating DWT with independent component analysis (ICA) was recently proposed and applied to structural damage detection. This method shows excellent performance in identifying sudden stiffness change for linear structures [6]. Independent component analysis (ICA) [30] is a computational technique for feature extraction in signal processing, which can separate a multivariate observed signal to statistically independent sources. It can be specially described as(15)wXt=AwSt,where wSt=wS1t,…,wSmt represents an n-dimensional independent source matrix and A represents an m × n (m ≥ n) mixing matrix.The principle of ICA can be summarized as an optimization process to search for proper estimation of the inverse ofA (i.e., demixing matrix W) such that the source signal matrix wS can be recovered by(16)wSt=WwXt.Therefore, there are two key issues in the ICA: (1) the proper objective function in order to decide whether the obtained source signals are statistically independent; (2) the effective algorithm to implement the optimization of the objective function.The Fast ICA proposed by Hyvärinen and Oja [31] is one of the most efficient and popular algorithms in revised ICA. It is based on a fixed-point iteration scheme maximizing non-Gaussianity as the objective function.To extract the damage information hidden in the noisy wavelet-domain signals, the detailed signal matrixwX is processed by using the Fast ICA, thereby nIC independent components ICj (j = 1, 2, …, nIC) can be obtained. Subsequently, the feature independent component (FIC) which contains the spike property can be extracted, and the position of the spike in the time history indicates the damage occurrence instant.


### 4.1.3. Initial Localization of Damage

After conducting the Fast ICA, the wavelet detailed signal matrixwX can be reconstructed by(17)wX=A×IC=∑i=1naj⋅ICj,where aj = [a1j, …, aij, …, amj]T is the mixing vector which is called as the source distribution vector (SDV). Herein, the element aij that is defined as the source distribution factor (SDF) represents the proportion of ICj distributed in wXi; ICj is the j-th independent component extracted by the Fast ICA.When ICj contained structural singularity information, it is called the feature independent component (FICk k = 1, 2, …, h). If aij is the largest among its SDV, it indicates its corresponding location, where the response measured in the i-th sensor location Xi contains the most FIC, is the possibly damaged location.The SDV can be normalized further by equation (18). Thus, the structural damage can be initially localized by the maximum among the NSDV of the FIC.(18)NSDF=NSDF1,NSDF2,⋯,NSDFi,NSDFmT,NSDFi=aij∑i=1maij,i=1,2,…,m.


## 4.1.1. Separation of Wavelet Detailed Components

To complete the damage novelty detection, the time-series dynamic responses is preprocessed by DWT in order to separate the wavelet detailed components. DWT [29] is a useful tool for time-frequency analysis, which achieves multiresolution analysis of a signal by decomposing it into high-frequency (detail) and low-frequency (approximation) components at each level. If a signal is decomposed into LL levels, it can be reconstructed by(14)X¨t=∑ll=1LLX¨dllt+X¨aLLt=∑ll=1LL∑kk=−∞∞dkk,llψkk,llt+∑kk=−∞∞ckk,LLϕkk,LLt,where X¨dllt is the detailed component at the ll-th scale level; X¨aLLt is the approximated component at the LL-th scale level; kk and ll are the translation and scale parameters, respectively; dkk,ll and ψkk,llt are the detail coefficient and wavelet function at the ll-th scale level, respectively; and ckl,LL and ϕkk,LLt are the approximation coefficient and scaling function at the LL-th scale level, respectively.Two parameters should be predetermined using DWT to detect the singularity (or discontinuity points) of a signal. The first one is selecting the proper wavelet basis and the other one is the decomposition level involving in DWT. After determining the parameters, DWT is, respectively, performed in the structural dynamic responsesXi (i = 1, 2, …, m; m is the total number of sensors) measured from each sensor location. Consequently, the wavelet detailed components wXi on a certain scale which contain the most damage information can be separated, respectively. So, the wavelet detailed signal matrix can be built as wX = [wX1, wX2 , …, wXm]T.


## 4.1.2. Extraction of Feature Independent Component

A new singularity analysis method called DWT-ICA by integrating DWT with independent component analysis (ICA) was recently proposed and applied to structural damage detection. This method shows excellent performance in identifying sudden stiffness change for linear structures [6]. Independent component analysis (ICA) [30] is a computational technique for feature extraction in signal processing, which can separate a multivariate observed signal to statistically independent sources. It can be specially described as(15)wXt=AwSt,where wSt=wS1t,…,wSmt represents an n-dimensional independent source matrix and A represents an m × n (m ≥ n) mixing matrix.The principle of ICA can be summarized as an optimization process to search for proper estimation of the inverse ofA (i.e., demixing matrix W) such that the source signal matrix wS can be recovered by(16)wSt=WwXt.Therefore, there are two key issues in the ICA: (1) the proper objective function in order to decide whether the obtained source signals are statistically independent; (2) the effective algorithm to implement the optimization of the objective function.The Fast ICA proposed by Hyvärinen and Oja [31] is one of the most efficient and popular algorithms in revised ICA. It is based on a fixed-point iteration scheme maximizing non-Gaussianity as the objective function.To extract the damage information hidden in the noisy wavelet-domain signals, the detailed signal matrixwX is processed by using the Fast ICA, thereby nIC independent components ICj (j = 1, 2, …, nIC) can be obtained. Subsequently, the feature independent component (FIC) which contains the spike property can be extracted, and the position of the spike in the time history indicates the damage occurrence instant.


## 4.1.3. Initial Localization of Damage

After conducting the Fast ICA, the wavelet detailed signal matrixwX can be reconstructed by(17)wX=A×IC=∑i=1naj⋅ICj,where aj = [a1j, …, aij, …, amj]T is the mixing vector which is called as the source distribution vector (SDV). Herein, the element aij that is defined as the source distribution factor (SDF) represents the proportion of ICj distributed in wXi; ICj is the j-th independent component extracted by the Fast ICA.When ICj contained structural singularity information, it is called the feature independent component (FICk k = 1, 2, …, h). If aij is the largest among its SDV, it indicates its corresponding location, where the response measured in the i-th sensor location Xi contains the most FIC, is the possibly damaged location.The SDV can be normalized further by equation (18). Thus, the structural damage can be initially localized by the maximum among the NSDV of the FIC.(18)NSDF=NSDF1,NSDF2,⋯,NSDFi,NSDFmT,NSDFi=aij∑i=1maij,i=1,2,…,m.


## 4.2. Refined Damage Assessment

First, the original time-series responses are divided intoks + 1 segments as X¨t1,X¨t2,…,X¨tks+1 according to the damage instant identified through FICks, ks = 1, 2, ..., h.Second, IMPSCO is performed on theks + 1 segment X¨t1,X¨t2,…,X¨tks+1, respectively; thus, the damage location and extent can be identified precisely according to the stiffness variation at each time period.


## 5. Numerical Study

### 5.1. Model Description

In order to validate the applicability and effectiveness of the proposed two-stage damage detection strategy for nonlinear structure, the nonlinear model presented in Section 2 was subjected to El-Centro earthquake excitation and Kobe earthquake excitation, respectively, at the base, as shown in Figure3. The original acceleration responses are presented in Figures 10 and 11. As shown in Table 9, two damage cases with different damage novelty, location, and extent were simulated and discussed here. To approach the operating environment, Gaussian white noise with SNR=30dB is added to all the responses and the earthquake excitations as equation (10).Figure 10
Original acceleration responses of the simulated structure under El-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 11
Original acceleration responses of the simulated structure under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Table 9
Damage scenarios.
Damage scenario
Damage novelty, location, and extent
Case 1
t = 10 s, k2 = 102.64 kN/m ⟶ 61.58 kN/m
Case 2
t = 10 s, k1 = 102.64 kN/m ⟶ 71.85 kN/m
t = 20 s, k2 = 102.64 kN/m ⟶ 71.85 kN/m


### 5.2. Damage Detection Model

#### 5.2.1. The First Stage: Detect Structural Response Abnormality and Initially Locate Damage Location

Just as Section4.1, the measured acceleration responses Ai (i = 1, 2, 3) were decomposed in single level by wavelet basis function db4, respectively. The corresponding detailed components wXi (i = 1, 2) were then selected and built as the wavelet detailed signal matrix wX = [wX1, wX2 ]T. The Fast ICA is performed on wX to obtain the two independent components (ICjj=1,2) and FICs. According to the position of spike in the FIC, the structural response novelty was identified, and the corresponding NSDV of the FIC was calculated via equation (18), and thus the damage location was initially estimated.


#### 5.2.2. The Second Stage: Time-Varying Parameter Identification

The original acceleration responsesAi (i = 1, 2) were divided into several segments according to the obtained FICs, and then IMPSCO was employed to identify the structural parameters in each segment, respectively. In this example, the structural mass was assumed as known in advance; thus, the parameter vector to be optimized was θ=k1,k2,c1,c2,n,γ,α,β,Y. The search range of θ was within the interval of [0.5b, 2b], and b is the theoretical value of θ. Additionally, the IMPSCO parameters were set as listed in Table 2. Meanwhile, because the IMPSCO algorithm is a probabilistic optimization algorithm, the identification process was performed for ten times and the average values were regarded as final results to eliminate the influence of randomness. Finally, the structural damages were localized and quantified accurately depending on the stiffness variation in each time period.


### 5.3. Results and Discussion

#### 5.3.1. Primary Damage Assessment

The results of wavelet detailed component separation are shown in Figures12 and 13, and the feature independent components (FICs) extracted are shown in Figures 14 and 15, respectively.Figure 12
Wavelet detailed components of the simulated structure under E1-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 13
Wavelet detailed components of the simulated structure under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 14
Independent component of the simulated structure under El-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 15
Independent component of the simulated structure under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)It is seen that significant damage information can hardly be directly identified from the wavelet detailed components when the structure is subjected to the earthquake excitation. After the Fast ICA processing as shown in Figures14 and 15, the FIC containing the spike is extracted obviously. More specifically, IC2 is the FIC in case 1, and it indicates that the structure is probably damaged at t = 10 s because the spike occurs at t = 10 s under both earthquake motions. Similarly, in case 2, IC2 can be used to estimate the occurring time of damage. It is found from Figure 15, as a consequence, that the structure is probably damaged at t = 10 s and t = 20 s, respectively.After detecting structural response novelty, the damage can then be initially localized by the NSDV. For single-damage pattern, case 1, we utilized NSDV to initially locate structural damage in case 1 under both earthquake excitations. The results are shown in Figure16. It can be seen that the maximum of NSDF corresponding to IC2 occurs at the second story when the structure was subject to El-Centro and Kobe earthquake excitation, respectively. It indicates that the second story is probably damaged location when t = 10 s. It can be seen from Figure 15 (b) that for mutidamage pattern, Case 2, damage abnormal information contained in IC2. Therefore, it is difficult to locate structural damage. All in all, structural response novelty can be identified correctly, and for single damage pattern, damage can be located initially.Figure 16
NSDV of the FIC under different earthquake waves. (a) EI-Centro earthquake. (b) Kobe earthquake.
(a)
(b)


#### 5.3.2. Refined Damage Assessment

When the primary damage assessment was completed, the IMPSCO was then used to implement segmentation parameters optimization and refined damage assessment. For case 1, the identification process was divided into two segments, namely, 0∼10 s and 10∼30 s; correspondingly, for case 2, the identification process is divided into three segments, namely, 0∼10 s, 10∼20 s, and 20∼30 s. Furthermore, the IMPSCP is used to identify the structural stiffnesses and coefficients as well as the parameters in the Bouc–Wen model, such as stiffness (A), loop fatness (β), loop pinching (γ), and abruptness parameters (n), and the results are shown in Figures 17 and 18. Tables 10 and 11 list the errors between identification results and theoretical values under different earthquake waves. It is found that the identification accuracy is very high for structural stiffnesses. More specifically, the maximum errors for the structural stiffness are 0.93% and 1.08% under EI-Centro earthquake wave, respectively. By contrast, the maximum errors are 0.51% and 1.25% under Kobe earthquake wave, respectively. On the contrary, except story stiffness, the rest of identified parameters have lower accuracy, but the maximum errors are no more than 12% and 14%, under two different earthquake waves. These error values basically met the requirements of engineering needs. In addition, difference in the maximum error values under the different ground motions used is due to differences in the structural responses induced by different earthquake damage.Figure 17
Structural parameters’ identification results when structure is subject to EI-Centro earthquake. (a) Case 1. (b) Case 2.
(a)
(b)Figure 18
Identification results under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Table 10
Relative errors between identification results and theoretical values under El-Centro earthquake wave.
Case
Time (s)
Relative error (%)
k
1
k
2
ζ
1
ζ
2
Y
α
β
γ
n
Case 1
0∼10
0.93
0.12
6.50
8.50
0.29
10.00
6.00
8.00
4.00
10∼30
0.89
0.10
6.50
2.00
4.00
10.00
0.00
2.00
2.00
Case 2
0∼10
0.49
0.35
8.00
8.00
5.07
10.00
8.00
8.00
4.00
10∼20
0.96
0.23
6.00
9.00
0.02
10.00
8.00
12.00
4.50
20∼30
1.08
0.07
3.00
8.00
5.49
10.00
8.00
8.00
5.50Table 11
Relative errors between identification results and theoretical values under Kobe earthquake wave.
Case
Time (s)
Relative error (%)
k
1
k
2
ζ
1
ζ
2
Y
α
β
γ
n
Case 1
0∼10
0.15
0.10
7.00
8.50
6.82
0.00
12.00
12.00
2.00
10∼30
0.51
0.29
6.00
7.50
4.91
0.00
14.00
12.00
3.00
Case 2
0∼10
1.25
0.14
11.00
8.00
5.93
0.00
6.00
14.00
2.50
10∼20
1.10
0.29
6.00
9.00
6.91
0.00
12.00
14.00
1.50
20∼30
1.17
1.04
6.00
7.50
4.60
0.00
12.00
10.00
2.00All in all, it can be seen that the identification results of the structural stiffnesses agree well with the theoretical values in both cases. This indicates that the precise damage localization and accurate assessment of damage extent can be realized in the second stage. Despite the rest of identified parameters with certain errors, these error values are good within modeling and construction errors. This implies that the proposed two-stage identification method is applicable to detecting sudden stiffness change in the acceleration responses for the nonlinear structures.


#### 5.3.3. Discussion

For the sake of further evaluating the proposed method, identified and actual hysteresis loops for the structure under different ground motions are shown in Figures19 and 20, respectively. In addition, identified and actual velocities and displacements for the structure subjected to different earthquakes waves are depicted in Figures 21 and 22, respectively. As shown in these figures, even for the worst damage pattern, case 2, with sudden changes in stiffness, the identified hysteresis loops, velocity, and displacements are in very good agreement with the actual values.Figure 19
Comparison of hysteresis loops under EI-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 20
Comparison of hysteresis loops under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 21
Comparison of velocities and displacements under El-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 22
Comparison of velocities and displacements under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)The identification results presented in this study are functionally equivalent to real-time results. Specifically, using software MATLAB R2015b, it took 8.96 s for each piecewise step of 10 s (500 Hz sampling rate) on a 3.20 GHz Intel® dual-core desktop machine to identify two stiffness values, damping, and the combined Bouc–Wen model parameters. In addition, the identification process at each time step only relies on the prior time step values. Hence, the proposed algorithm can be readily used as an online SHM method. In other words, the proposed strategy is a good choice for identifying sudden stiffness change in the acceleration response of a nonlinear hysteretic structure.


## 5.1. Model Description

In order to validate the applicability and effectiveness of the proposed two-stage damage detection strategy for nonlinear structure, the nonlinear model presented in Section 2 was subjected to El-Centro earthquake excitation and Kobe earthquake excitation, respectively, at the base, as shown in Figure3. The original acceleration responses are presented in Figures 10 and 11. As shown in Table 9, two damage cases with different damage novelty, location, and extent were simulated and discussed here. To approach the operating environment, Gaussian white noise with SNR=30dB is added to all the responses and the earthquake excitations as equation (10).Figure 10
Original acceleration responses of the simulated structure under El-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 11
Original acceleration responses of the simulated structure under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Table 9
Damage scenarios.
Damage scenario
Damage novelty, location, and extent
Case 1
t = 10 s, k2 = 102.64 kN/m ⟶ 61.58 kN/m
Case 2
t = 10 s, k1 = 102.64 kN/m ⟶ 71.85 kN/m
t = 20 s, k2 = 102.64 kN/m ⟶ 71.85 kN/m


## 5.2. Damage Detection Model

### 5.2.1. The First Stage: Detect Structural Response Abnormality and Initially Locate Damage Location

Just as Section4.1, the measured acceleration responses Ai (i = 1, 2, 3) were decomposed in single level by wavelet basis function db4, respectively. The corresponding detailed components wXi (i = 1, 2) were then selected and built as the wavelet detailed signal matrix wX = [wX1, wX2 ]T. The Fast ICA is performed on wX to obtain the two independent components (ICjj=1,2) and FICs. According to the position of spike in the FIC, the structural response novelty was identified, and the corresponding NSDV of the FIC was calculated via equation (18), and thus the damage location was initially estimated.


### 5.2.2. The Second Stage: Time-Varying Parameter Identification

The original acceleration responsesAi (i = 1, 2) were divided into several segments according to the obtained FICs, and then IMPSCO was employed to identify the structural parameters in each segment, respectively. In this example, the structural mass was assumed as known in advance; thus, the parameter vector to be optimized was θ=k1,k2,c1,c2,n,γ,α,β,Y. The search range of θ was within the interval of [0.5b, 2b], and b is the theoretical value of θ. Additionally, the IMPSCO parameters were set as listed in Table 2. Meanwhile, because the IMPSCO algorithm is a probabilistic optimization algorithm, the identification process was performed for ten times and the average values were regarded as final results to eliminate the influence of randomness. Finally, the structural damages were localized and quantified accurately depending on the stiffness variation in each time period.


## 5.2.1. The First Stage: Detect Structural Response Abnormality and Initially Locate Damage Location

Just as Section4.1, the measured acceleration responses Ai (i = 1, 2, 3) were decomposed in single level by wavelet basis function db4, respectively. The corresponding detailed components wXi (i = 1, 2) were then selected and built as the wavelet detailed signal matrix wX = [wX1, wX2 ]T. The Fast ICA is performed on wX to obtain the two independent components (ICjj=1,2) and FICs. According to the position of spike in the FIC, the structural response novelty was identified, and the corresponding NSDV of the FIC was calculated via equation (18), and thus the damage location was initially estimated.


## 5.2.2. The Second Stage: Time-Varying Parameter Identification

The original acceleration responsesAi (i = 1, 2) were divided into several segments according to the obtained FICs, and then IMPSCO was employed to identify the structural parameters in each segment, respectively. In this example, the structural mass was assumed as known in advance; thus, the parameter vector to be optimized was θ=k1,k2,c1,c2,n,γ,α,β,Y. The search range of θ was within the interval of [0.5b, 2b], and b is the theoretical value of θ. Additionally, the IMPSCO parameters were set as listed in Table 2. Meanwhile, because the IMPSCO algorithm is a probabilistic optimization algorithm, the identification process was performed for ten times and the average values were regarded as final results to eliminate the influence of randomness. Finally, the structural damages were localized and quantified accurately depending on the stiffness variation in each time period.


## 5.3. Results and Discussion

### 5.3.1. Primary Damage Assessment

The results of wavelet detailed component separation are shown in Figures12 and 13, and the feature independent components (FICs) extracted are shown in Figures 14 and 15, respectively.Figure 12
Wavelet detailed components of the simulated structure under E1-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 13
Wavelet detailed components of the simulated structure under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 14
Independent component of the simulated structure under El-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 15
Independent component of the simulated structure under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)It is seen that significant damage information can hardly be directly identified from the wavelet detailed components when the structure is subjected to the earthquake excitation. After the Fast ICA processing as shown in Figures14 and 15, the FIC containing the spike is extracted obviously. More specifically, IC2 is the FIC in case 1, and it indicates that the structure is probably damaged at t = 10 s because the spike occurs at t = 10 s under both earthquake motions. Similarly, in case 2, IC2 can be used to estimate the occurring time of damage. It is found from Figure 15, as a consequence, that the structure is probably damaged at t = 10 s and t = 20 s, respectively.After detecting structural response novelty, the damage can then be initially localized by the NSDV. For single-damage pattern, case 1, we utilized NSDV to initially locate structural damage in case 1 under both earthquake excitations. The results are shown in Figure16. It can be seen that the maximum of NSDF corresponding to IC2 occurs at the second story when the structure was subject to El-Centro and Kobe earthquake excitation, respectively. It indicates that the second story is probably damaged location when t = 10 s. It can be seen from Figure 15 (b) that for mutidamage pattern, Case 2, damage abnormal information contained in IC2. Therefore, it is difficult to locate structural damage. All in all, structural response novelty can be identified correctly, and for single damage pattern, damage can be located initially.Figure 16
NSDV of the FIC under different earthquake waves. (a) EI-Centro earthquake. (b) Kobe earthquake.
(a)
(b)


### 5.3.2. Refined Damage Assessment

When the primary damage assessment was completed, the IMPSCO was then used to implement segmentation parameters optimization and refined damage assessment. For case 1, the identification process was divided into two segments, namely, 0∼10 s and 10∼30 s; correspondingly, for case 2, the identification process is divided into three segments, namely, 0∼10 s, 10∼20 s, and 20∼30 s. Furthermore, the IMPSCP is used to identify the structural stiffnesses and coefficients as well as the parameters in the Bouc–Wen model, such as stiffness (A), loop fatness (β), loop pinching (γ), and abruptness parameters (n), and the results are shown in Figures 17 and 18. Tables 10 and 11 list the errors between identification results and theoretical values under different earthquake waves. It is found that the identification accuracy is very high for structural stiffnesses. More specifically, the maximum errors for the structural stiffness are 0.93% and 1.08% under EI-Centro earthquake wave, respectively. By contrast, the maximum errors are 0.51% and 1.25% under Kobe earthquake wave, respectively. On the contrary, except story stiffness, the rest of identified parameters have lower accuracy, but the maximum errors are no more than 12% and 14%, under two different earthquake waves. These error values basically met the requirements of engineering needs. In addition, difference in the maximum error values under the different ground motions used is due to differences in the structural responses induced by different earthquake damage.Figure 17
Structural parameters’ identification results when structure is subject to EI-Centro earthquake. (a) Case 1. (b) Case 2.
(a)
(b)Figure 18
Identification results under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Table 10
Relative errors between identification results and theoretical values under El-Centro earthquake wave.
Case
Time (s)
Relative error (%)
k
1
k
2
ζ
1
ζ
2
Y
α
β
γ
n
Case 1
0∼10
0.93
0.12
6.50
8.50
0.29
10.00
6.00
8.00
4.00
10∼30
0.89
0.10
6.50
2.00
4.00
10.00
0.00
2.00
2.00
Case 2
0∼10
0.49
0.35
8.00
8.00
5.07
10.00
8.00
8.00
4.00
10∼20
0.96
0.23
6.00
9.00
0.02
10.00
8.00
12.00
4.50
20∼30
1.08
0.07
3.00
8.00
5.49
10.00
8.00
8.00
5.50Table 11
Relative errors between identification results and theoretical values under Kobe earthquake wave.
Case
Time (s)
Relative error (%)
k
1
k
2
ζ
1
ζ
2
Y
α
β
γ
n
Case 1
0∼10
0.15
0.10
7.00
8.50
6.82
0.00
12.00
12.00
2.00
10∼30
0.51
0.29
6.00
7.50
4.91
0.00
14.00
12.00
3.00
Case 2
0∼10
1.25
0.14
11.00
8.00
5.93
0.00
6.00
14.00
2.50
10∼20
1.10
0.29
6.00
9.00
6.91
0.00
12.00
14.00
1.50
20∼30
1.17
1.04
6.00
7.50
4.60
0.00
12.00
10.00
2.00All in all, it can be seen that the identification results of the structural stiffnesses agree well with the theoretical values in both cases. This indicates that the precise damage localization and accurate assessment of damage extent can be realized in the second stage. Despite the rest of identified parameters with certain errors, these error values are good within modeling and construction errors. This implies that the proposed two-stage identification method is applicable to detecting sudden stiffness change in the acceleration responses for the nonlinear structures.


### 5.3.3. Discussion

For the sake of further evaluating the proposed method, identified and actual hysteresis loops for the structure under different ground motions are shown in Figures19 and 20, respectively. In addition, identified and actual velocities and displacements for the structure subjected to different earthquakes waves are depicted in Figures 21 and 22, respectively. As shown in these figures, even for the worst damage pattern, case 2, with sudden changes in stiffness, the identified hysteresis loops, velocity, and displacements are in very good agreement with the actual values.Figure 19
Comparison of hysteresis loops under EI-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 20
Comparison of hysteresis loops under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 21
Comparison of velocities and displacements under El-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 22
Comparison of velocities and displacements under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)The identification results presented in this study are functionally equivalent to real-time results. Specifically, using software MATLAB R2015b, it took 8.96 s for each piecewise step of 10 s (500 Hz sampling rate) on a 3.20 GHz Intel® dual-core desktop machine to identify two stiffness values, damping, and the combined Bouc–Wen model parameters. In addition, the identification process at each time step only relies on the prior time step values. Hence, the proposed algorithm can be readily used as an online SHM method. In other words, the proposed strategy is a good choice for identifying sudden stiffness change in the acceleration response of a nonlinear hysteretic structure.


## 5.3.1. Primary Damage Assessment

The results of wavelet detailed component separation are shown in Figures12 and 13, and the feature independent components (FICs) extracted are shown in Figures 14 and 15, respectively.Figure 12
Wavelet detailed components of the simulated structure under E1-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 13
Wavelet detailed components of the simulated structure under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 14
Independent component of the simulated structure under El-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 15
Independent component of the simulated structure under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)It is seen that significant damage information can hardly be directly identified from the wavelet detailed components when the structure is subjected to the earthquake excitation. After the Fast ICA processing as shown in Figures14 and 15, the FIC containing the spike is extracted obviously. More specifically, IC2 is the FIC in case 1, and it indicates that the structure is probably damaged at t = 10 s because the spike occurs at t = 10 s under both earthquake motions. Similarly, in case 2, IC2 can be used to estimate the occurring time of damage. It is found from Figure 15, as a consequence, that the structure is probably damaged at t = 10 s and t = 20 s, respectively.After detecting structural response novelty, the damage can then be initially localized by the NSDV. For single-damage pattern, case 1, we utilized NSDV to initially locate structural damage in case 1 under both earthquake excitations. The results are shown in Figure16. It can be seen that the maximum of NSDF corresponding to IC2 occurs at the second story when the structure was subject to El-Centro and Kobe earthquake excitation, respectively. It indicates that the second story is probably damaged location when t = 10 s. It can be seen from Figure 15 (b) that for mutidamage pattern, Case 2, damage abnormal information contained in IC2. Therefore, it is difficult to locate structural damage. All in all, structural response novelty can be identified correctly, and for single damage pattern, damage can be located initially.Figure 16
NSDV of the FIC under different earthquake waves. (a) EI-Centro earthquake. (b) Kobe earthquake.
(a)
(b)


## 5.3.2. Refined Damage Assessment

When the primary damage assessment was completed, the IMPSCO was then used to implement segmentation parameters optimization and refined damage assessment. For case 1, the identification process was divided into two segments, namely, 0∼10 s and 10∼30 s; correspondingly, for case 2, the identification process is divided into three segments, namely, 0∼10 s, 10∼20 s, and 20∼30 s. Furthermore, the IMPSCP is used to identify the structural stiffnesses and coefficients as well as the parameters in the Bouc–Wen model, such as stiffness (A), loop fatness (β), loop pinching (γ), and abruptness parameters (n), and the results are shown in Figures 17 and 18. Tables 10 and 11 list the errors between identification results and theoretical values under different earthquake waves. It is found that the identification accuracy is very high for structural stiffnesses. More specifically, the maximum errors for the structural stiffness are 0.93% and 1.08% under EI-Centro earthquake wave, respectively. By contrast, the maximum errors are 0.51% and 1.25% under Kobe earthquake wave, respectively. On the contrary, except story stiffness, the rest of identified parameters have lower accuracy, but the maximum errors are no more than 12% and 14%, under two different earthquake waves. These error values basically met the requirements of engineering needs. In addition, difference in the maximum error values under the different ground motions used is due to differences in the structural responses induced by different earthquake damage.Figure 17
Structural parameters’ identification results when structure is subject to EI-Centro earthquake. (a) Case 1. (b) Case 2.
(a)
(b)Figure 18
Identification results under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Table 10
Relative errors between identification results and theoretical values under El-Centro earthquake wave.
Case
Time (s)
Relative error (%)
k
1
k
2
ζ
1
ζ
2
Y
α
β
γ
n
Case 1
0∼10
0.93
0.12
6.50
8.50
0.29
10.00
6.00
8.00
4.00
10∼30
0.89
0.10
6.50
2.00
4.00
10.00
0.00
2.00
2.00
Case 2
0∼10
0.49
0.35
8.00
8.00
5.07
10.00
8.00
8.00
4.00
10∼20
0.96
0.23
6.00
9.00
0.02
10.00
8.00
12.00
4.50
20∼30
1.08
0.07
3.00
8.00
5.49
10.00
8.00
8.00
5.50Table 11
Relative errors between identification results and theoretical values under Kobe earthquake wave.
Case
Time (s)
Relative error (%)
k
1
k
2
ζ
1
ζ
2
Y
α
β
γ
n
Case 1
0∼10
0.15
0.10
7.00
8.50
6.82
0.00
12.00
12.00
2.00
10∼30
0.51
0.29
6.00
7.50
4.91
0.00
14.00
12.00
3.00
Case 2
0∼10
1.25
0.14
11.00
8.00
5.93
0.00
6.00
14.00
2.50
10∼20
1.10
0.29
6.00
9.00
6.91
0.00
12.00
14.00
1.50
20∼30
1.17
1.04
6.00
7.50
4.60
0.00
12.00
10.00
2.00All in all, it can be seen that the identification results of the structural stiffnesses agree well with the theoretical values in both cases. This indicates that the precise damage localization and accurate assessment of damage extent can be realized in the second stage. Despite the rest of identified parameters with certain errors, these error values are good within modeling and construction errors. This implies that the proposed two-stage identification method is applicable to detecting sudden stiffness change in the acceleration responses for the nonlinear structures.


## 5.3.3. Discussion

For the sake of further evaluating the proposed method, identified and actual hysteresis loops for the structure under different ground motions are shown in Figures19 and 20, respectively. In addition, identified and actual velocities and displacements for the structure subjected to different earthquakes waves are depicted in Figures 21 and 22, respectively. As shown in these figures, even for the worst damage pattern, case 2, with sudden changes in stiffness, the identified hysteresis loops, velocity, and displacements are in very good agreement with the actual values.Figure 19
Comparison of hysteresis loops under EI-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 20
Comparison of hysteresis loops under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 21
Comparison of velocities and displacements under El-Centro earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)Figure 22
Comparison of velocities and displacements under Kobe earthquake wave. (a) Case 1. (b) Case 2.
(a)
(b)The identification results presented in this study are functionally equivalent to real-time results. Specifically, using software MATLAB R2015b, it took 8.96 s for each piecewise step of 10 s (500 Hz sampling rate) on a 3.20 GHz Intel® dual-core desktop machine to identify two stiffness values, damping, and the combined Bouc–Wen model parameters. In addition, the identification process at each time step only relies on the prior time step values. Hence, the proposed algorithm can be readily used as an online SHM method. In other words, the proposed strategy is a good choice for identifying sudden stiffness change in the acceleration response of a nonlinear hysteretic structure.


## 6. Conclusions and Remarks

This paper proposes a two-stage damage detection strategy by combining the DWT-FastICA technology and the IMPSCO algorithm for detecting sudden stiffness change in the accelerations of nonlinear hysteretic structures. Numerical simulation and experiment are used to validate the efficiency of the proposed methods. Some conclusions can be drawn as follows:(1)
The proposed strategy can not only locate the structural damage initially and quickly but also quantify the damage severity accurately. It is noted that the proposed strategy is implemented only by using any kinds of structural time-series responses and does not require structural displacement measurements which are typically difficult to acquire or reasonably estimate.(2)
The IMPSCO algorithm is sui` to identify parameters in the modified Bouc–Wen model simulated by a MR damper, and the reconstruction force-displacement data can match well with relevant experimental data. The identification process at each time step based on IMPSCO is computationally efficient for sudden stiffness change in the acceleration response.(3)
The integration of DWT and Fast ICA can not only detect structural response novelty but also primary locate structural damage induced by stiffness change in time domain.(4)
The proposed strategy is applicable and effective for identifying key structural parameters including stiffness, damping, and the Bouc–Wen hysteretic model parameters with sudden stiffness change in the acceleration response. These parameters are directly related to well-recognized damage metrics.(5)
Proof-of-method simulations of a realistic nonlinear case-study structure, subjected to a suite of two different ground motions, show that the proposed method is well capable of identifying structural parameters to within 1.25% of the actual as-modeled values for stiffness.In summary, it is effective and rapid to identify sudden stiffness change in the acceleration response of the two-story nonlinear Bouc–Wen model, and it still needs more experiments and structures to validate its efficiency in the future.



---

*Source: 3824216-2020-01-31.xml*