# Whale Optimization Algorithm-Based LQG-Adaptive Neuro-Fuzzy Control for Seismic Vibration Mitigation with MR Dampers

Authors: Xiufang Lin; Weiqing Lin
Journal: Shock and Vibration
 (2022)

---

## Abstract

Design of effective control strategies to protect structural buildings from seismic hazards is gaining increasing attention. In this paper, an intelligent semiactive control strategy, which combines linear-quadratic-Gaussian (LQG), whale optimization algorithm (WOA), and adaptive neuro-fuzzy inference system (ANFIS) strategy is designed to mitigate structural vibration by using magnetorheological (MR) dampers, here known as WLQG-ANFIS control. Firstly, considering that the performance of the LQG control for the structures under seismic excitations strongly depends on skill and experience of the experts in determining the weighting matrix of the feedback gain, a WOA with constraints is adopted to optimally design the LQG control due to the strong ability to avoid local optima. Secondly, the inverse dynamic model of the MR damper is developed with the ANFIS technique which combines the fuzzy system and artificial neural network. Finally, according to the active force estimated by the proposed LQG control, the developed inverse model is utilized to calculate the MR damper command signal. Numerical analyses are conducted for a 10-story structural building installed with MR dampers under earthquake excitation. Through a test called analysis of variance in optimizing the LQG control, the superiority of WOA over other three meta-heuristic algorithms, i.e. GA, DE, and ABC, in terms of convergence performance and robustness is validated. The proposed WLQG-ANFIS control can reduce the maximum displacement, interstory drift and acceleration responses by 64%, 65%, and 53%, respectively, which can achieve better overall performance than on-off control, LQG-ANFIS control, fuzzy control, and H∞-ANFIS control. Moreover, this control strategy exhibits the ability to handle the uncertainties in structural parameters as well as seismic excitations.


---

## Body

## 1. Introduction

Nowadays, it has gaining more and more attention in civil engineering to mitigate the structural vibration under earthquake excitations which could endanger the safety of people’s lives and property. Generally, there are three kinds of structural control systems, i.e. passive, active, and semiactive control systems. The semiactive control system combines the advantages of the other two systems so that it is highly reliable and significantly adaptable, which makes this control system as a very promising one [1].The magnetorheological (MR) damper is a kind of intelligent semi-active control device which can generate large force with very low power and is fail-safe due to the excellent characteristics of the MR fluid. Therefore, the MR damper has attracted wide attention in the field of vibration control such as civil structures [2, 3], vehicle system [4], offshore platform [5], and washing machines [6].Until now, there are many kinds of semiactive control methods for the MR damper-based vibration control. For example, decentralized bang-bang control [7], Lyapunov method [8] and modulated homogeneous friction control [9] are traditional control methods, most of which supplied either zero or the maximum value for the command signal of MR dampers. H∞ control [10] and sliding mode control [11, 12] are robust control methods which have the capability to maintain stability and performance specifications despite the uncertainties of the system. Neural network control [13] and fuzzy logic control (FLC) [14, 15] are intelligent control algorithms which can handle high non-linearities and heuristic knowledge. Do et al. [16] proposed a new direct adaptive fuzzy controller for the damping force tracking control of MR damper system. They also proposed a robust control strategy using adaptive fuzzy neural network for vibration control systems with dead-zone band and time delay under severe disturbance [17]. An new hybrid optimal controller including several techniques such as fuzzy model, optimal control, sliding mode control and H-infinity was also designed and applied to a vehicle seat suspension for vibration control [18]. Besides, linear quadratic optimal control systems including linear quadratic regulator (LQR) control and linear quadratic Gaussian (LQG) control have been widely used. A modified linear matrix inequality-based LQR was proposed for the semiactive control with MR dampers [19]. A classical LQR controller was selected as the primary controller to calculate desired control force based on the structural response of the MR elastomer base isolator [20]. Lee et al. [21] developed an optimal path tracking control based on LQG for autonomous vehicles. However, a major drawback of the LQR control is that it requires full-state feedback which cannot be easily measured in the real application, especially when the controlled object is a high-rise structural building. Therefore, as an output feedback type, LQG control which can estimate full states through Kalman filter is more favorable to engineering application compared with LQR control. Wanf and Dyke [22] proposed a modal-based LQG to suppress the seismic responses for smart base isolation system with MR dampers.However, similar to the LQR control, performance of the LQG highly requires proper selecting the weighting matrix parameters of feedback gain, and yet the chosen of these parameters for the vibration control of the seismically excited structures has no general criteria to follow and relies on the experts’ skill. This will result in low efficiency in controller design and possible unsatisfactory control effects. Till now, it has become a trend to tune the gain matrix parameters of the linear quadratic optimal controllers with meta-heuristic algorithms (MAs). Mansoor et al. [23] employed genetic algorithm (GA) for optimally choosing the optimal weighting matrix parameters of the LQR control for hovercraft. The results showed that the step response is very satisfactory. To design the LQG autonomous controllers for surgical robots, a particle swarm optimization (PSO) was utilized to optimally select the weighting matrices of the controller. The results verified that this PSO-based LQG control yielded higher convergence speed and accuracy than that whose weighting matrices selected by trial and error [24]. However, as far as we know, the MA-optimized LQG controller has not been considered for the vibration control with MR dampers.As a new kind of meta-heuristic algorithm, whale optimization algorithm (WOA) was proposed by Mirjalili and Lewis at the first time to solve mathematical and engineering problems [25]. The results demonstrates that this algorithm is very competitive among several meta-heuristic methods, e.g. differential evolution (DE) and gravitational search algorithm. The working principle of the WOA imitates the social behavior of humpback whales. Besides, in comparison to GA, WOA has no evolution operations and easy to be conducted [26]. Therefore, WOA combines the advantages of the genetic-based memetic algorithm and the social behavior-based PSO. It should be noted that WOA conducts exploration and two kinds of exploitations, which prevents this algorithm from getting stuck in local optima and makes it converge fast.Recently, WOA has found its application in various engineering fields [27]. Specifically, Aljarah et al. adopted WOA to train multilayer perception neural networks and proved that as for the majority of datasets, WOA outperformed back-propagation algorithm, GA, artificial bee colony (ABC), and ant colony optimization on in terms of accuracy and speed of convergence [28]. WOA was also used to solve the reactive power planning problem [29]. The superiority of the WOA was verified through the comparison with other meta-heuristic techniques, such as DE and grey wolf optimization. Hasanien and Hany [30] applied WOA to design the proportional-integral controllers of the Photovoltaic system. Analysis results showed that the proposed controller outperformed the controller optimized with generalized reduced gradient algorithm. Lin and Lin [31] used WOA to simultaneously optimize the damper allocation and fuzzy controller for the vibration control of two connected adjacent buildings. The comparison results show that WOA performed better than GA and PSO in terms of local optima avoidance and robustness. What’s more, WOA was applied in many other fields of engineering optimization such as supply of the water reservoir demand [32], fuzzy logic controller design for the wind generators [33], structural damage identification [34], and chaotic map for enhanced image encryption [35]. In short, it has been validated that WOA is very promising and encouraging.On the other hand, in the MR damper-based semiactive control, the damping force cannot be directly regulated on account of the damper’s multi-physics and intrinsic nonlinearity. Therefore, generally it is necessary to develop an MR damper controller to determine the input signal for accurately controlling the damping force, so as to track the active force calculated by the system controller. There are mainly two kinds of methods to design the MR damper controllers. One method is using the force feedback control with force sensors, where voltage prediction approach, like clipped voltage law [36] is used. But with this method damping force cannot track the active force precisely because it cannot output continual voltage value to the MR damper. Besides, the complicacy and the cost of the control system will inevitably be increased due to the requirements of extraforce sensors.The other method to design the damper controller is developing an inverse dynamic model, including parametric and nonparametric methods. The parametric models are either not accurate enough due to the simplicity of some forward models or difficult to derive due to the high nonlinearity of some forward models, e.g., phenomenological model [37]. As an intelligent nonparametric modeling strategy, a new neural network is developed to characterize the inverse model for MR elastomer base isolator and the effectiveness was proved by experimental results [20]. Besides, adaptive neuro-fuzzy inference system (ANFIS) proposed by Jang [38] combines the fuzzy system and artificial neural network. It possesses excellent reasoning capability of the fuzzy system as well as learning skills and computational power of the artificial neural network. Recently, this modeling strategy has been proved to be very promising to develop the MR damper inverse model. Lin et al. [39] utilized the ANFIS technique to characterize an MR damper inverse model for the vibration control of a seismically-excited structural building. The results proved that by use of the ANFIS inverse model, the desired command voltage matched well with the predicted command voltage, and the seismic responses can be significantly mitigated.In this study, a semiactive control strategy using MR dampers integrating WOA-based LQG with ANFIS is proposed to improve the seismic performance of the structural building. To estimate the optimal active force of the control system, motivated by the obvious advantages of WOA, this optimization algorithm is utilized to design the optimal feedback gain of the LQG controller. Besides, based on the optimal active force, to generate precise damping force, the inverse model is developed by the ANFIS approach to predict the input signal of the MR damper. For investigating the performance of the proposed control strategy, a building installed with dampers undergoing earthquake excitations is considered and different control methods are compared in the simulation study. Besides, the control results with the consideration of uncertainties in structural parameters as well as seismic excitations are also analyzed. In summary, the major contributions of this paper can be highlighted as follows:(1)
A parameter identification approach for the LQG control using the WOA is proposed, which can effectively optimize the controller parameters. To the best of our knowledge, there has no literature published about this problem.(2)
The WOA-based LQG controller is incorporated with the MR damper inverse model developed using ANFIS for the semi-active control. As far as we know, this semiactive structural control strategy has not been studied before.(3)
An analysis of variance (ANOVA) test is employed for statistically evaluating the performance difference between the WOA and other three algorithms.(4)
The proposed WLQG-ANFIS semiactive control can decrease the structural responses more effectively than the other four semiactive control methods.The rest of this paper is organized as follows. Section2 presents the mathematical calculation model of a structural building installed with MR dampers under earthquake excitations, and then derives the corresponding motion equation. Section 3 provides an overview of the WOA, derives the weighting matrices of LQG based on the desired control performance and illustrates the design of the proposed WOA-based LQG control strategy. In Section 4, the ANFIS technology is utilized to develop the MR damper inverse model. The simulation test results are given in Section 5. The conclusion and future work are stated in Section 6.


## 2. Modeling the Controlled System

### 2.1. Motion Equation and State-Space of the System

Assume a steel-frame structure withm MR dampers having n stories excited by earthquake acceleration x¨gt, the motion equation can be shown as follows:(1)Mx¨t+Cx˙t+Kxt=Γft+MΛx¨gt,where the vector f(t) = [f1(t), f2(t),…, fm(t)]T represents the damping forces described in the next subsection; M,C,K∈Rn×n indicate mass, damping, and stiffness matrices of the structure, respectively. Γ∈Rn×m is the location of the m dampers; Λ denotes the influence of seismic excitation; and x(t) = [x1(t), x2(t),…, xn(t)]T is a displacement vector relative to ground.Ifzt=x˙tT,xtTT represents state variable, then we can rewrite (1) as(2)z˙t=Azt+B1ft+B2x¨gt,where(3)A=−M−1C−M−1KI0,B1=M−1Γ0,B2=Λ0,where I is the unit matrix; 0 is the zero matrix (vector).Ifx¨at denotes the absolute acceleration, then the system output denoted by yt=x¨atT,x˙tT,xtTT is expressed by(4)yt=Czzt+D1ft,where(5)Cz=−M−1C−M−1KI00I,D1=M−1Γ00.In order to conveniently conduct the simulation calculation, we describe (2) and (4) in the form of standard state-space as follows:(6)z˙t=Azt+But.(7)yt=Czzt+Dut,where B = [B1B2], D = [D10]; ut=ftT,x¨gtT is the system input.


### 2.2. Forward Dynamic Model of the MR Damper

The phenomenological model is used to portray the MR damper forward dynamics by which the damping force is calculated as follows [37].(8)f=c1y˙+k1x−x0,where(9)y˙=1c0+c1αz+c0x˙+k0x−y,where z is the evolutionary variable which is expressed as(10)z˙=−γx˙−y˙zzn−1−βx˙−y˙zn+Ax˙−y˙.In Equations (8)–(10), x is the displacement of the damper, and ẋ is the velocity of the damper. C0 is viscous damping at the high velocity. C1 is viscous damping at the low velocity. K0 denotes the stiffness constant at the high velocity. K1 presents the accumulator stiffness. X0 is the initial displacement for the spring element. α denotes a factor related to the hysteresis loop. Yielding element coefficients to adjust the hysteresis loop are represented by β, A, and γ. The following three equations are used to calculate α, c0, and c1, respectively.(11)α=αu=αa+αbu,(12)c0=c0u=c0a+c0bu,(13)c1=c1u=c1a+c1bu,where u is the function of the command voltage v shown as bellow(14)u˙=−ηu−v,where η is a constant.We adopt the MR damper which can generate the maximum force 1000 kN with the maximum command voltage 9 V in this study. All the parameters of the phenomenological model in equations (8)–(14) for this MR damper are listed in [40].


## 2.1. Motion Equation and State-Space of the System

Assume a steel-frame structure withm MR dampers having n stories excited by earthquake acceleration x¨gt, the motion equation can be shown as follows:(1)Mx¨t+Cx˙t+Kxt=Γft+MΛx¨gt,where the vector f(t) = [f1(t), f2(t),…, fm(t)]T represents the damping forces described in the next subsection; M,C,K∈Rn×n indicate mass, damping, and stiffness matrices of the structure, respectively. Γ∈Rn×m is the location of the m dampers; Λ denotes the influence of seismic excitation; and x(t) = [x1(t), x2(t),…, xn(t)]T is a displacement vector relative to ground.Ifzt=x˙tT,xtTT represents state variable, then we can rewrite (1) as(2)z˙t=Azt+B1ft+B2x¨gt,where(3)A=−M−1C−M−1KI0,B1=M−1Γ0,B2=Λ0,where I is the unit matrix; 0 is the zero matrix (vector).Ifx¨at denotes the absolute acceleration, then the system output denoted by yt=x¨atT,x˙tT,xtTT is expressed by(4)yt=Czzt+D1ft,where(5)Cz=−M−1C−M−1KI00I,D1=M−1Γ00.In order to conveniently conduct the simulation calculation, we describe (2) and (4) in the form of standard state-space as follows:(6)z˙t=Azt+But.(7)yt=Czzt+Dut,where B = [B1B2], D = [D10]; ut=ftT,x¨gtT is the system input.


## 2.2. Forward Dynamic Model of the MR Damper

The phenomenological model is used to portray the MR damper forward dynamics by which the damping force is calculated as follows [37].(8)f=c1y˙+k1x−x0,where(9)y˙=1c0+c1αz+c0x˙+k0x−y,where z is the evolutionary variable which is expressed as(10)z˙=−γx˙−y˙zzn−1−βx˙−y˙zn+Ax˙−y˙.In Equations (8)–(10), x is the displacement of the damper, and ẋ is the velocity of the damper. C0 is viscous damping at the high velocity. C1 is viscous damping at the low velocity. K0 denotes the stiffness constant at the high velocity. K1 presents the accumulator stiffness. X0 is the initial displacement for the spring element. α denotes a factor related to the hysteresis loop. Yielding element coefficients to adjust the hysteresis loop are represented by β, A, and γ. The following three equations are used to calculate α, c0, and c1, respectively.(11)α=αu=αa+αbu,(12)c0=c0u=c0a+c0bu,(13)c1=c1u=c1a+c1bu,where u is the function of the command voltage v shown as bellow(14)u˙=−ηu−v,where η is a constant.We adopt the MR damper which can generate the maximum force 1000 kN with the maximum command voltage 9 V in this study. All the parameters of the phenomenological model in equations (8)–(14) for this MR damper are listed in [40].


## 3. WOA-Based LQG Controller Design for Vibration Control

### 3.1. Design of LQG Control Strategy

In this study, the LQG controller is designed to calculate the desired controlling force. It consists of two parts, i.e. regulator and observer filter which can calculate feedback force and estimate the full states separately.


#### 3.1.1. Regulator Design

For the sake of security of the structure and the occupants’ comfort, the structural responses including maximum displacement and acceleration are required to be reduced. Besides, the desired controlling force should be limited in a reasonable range. Therefore, in this study, the index of control performance of the LQG is defined as(15)J=limx⟶∞1T∫0Tq1xmax2+q2x¨max2+∑i=1mrifdidt,where q1 and q2 are the weighting coefficients of the maximum relative displacement xmax and maximum absolute acceleration x¨max, respectively. For the sake of convenience, the relative displacement and absolute acceleration are called displacement and acceleration in short, respectively, in the subsequent discussion. The resultant feedback force fd which is a diagonal matrix consisting of fdi, also called the active force in the semi-active control, is used as one of the inputs for ANFIS inverse model to predict the command voltage (see Section 4.2).ri is the weighting coefficient of the feedback force fdi. M is the total number of the feedback forces, which is equal to the number of MR dampers. It should be noted that q1, q2, and ri are the parameters to be determined by WOA.The core of the control is to design a regulator so as to minimize the performance indexJ, which is rewritten as(16)J=limx⟶∞1T∫0TzTQz+fdTRfddt,where Q is the weighting matrix of state variable vector z defined in (2); R is the weighting matrix of the resultant feedback force fd defined in (15).According to the optimal control law, the resultant controlling force is written as(17)fdt=−12R−1BTPZ=−GZt,where G is the feedback gain, and P can be calculated by the following Riccati equation:(18)-PA-ATP+PBR−1BTP−Q=0,where A and B are the state equation matrices defined in (2) and (6), respectively, both of which are determined by the structural parameters.Finally, to calculate the optimal feedback gainG, the key step is to determine its weighting matrices Q and R. Based on (15) and (16), these weighting matrices are deduced as(19)Q=0000000000000000000⋱0000000q2Cn,n−1Mn,n2q2Cn,n−1Cn,nMn,n2000q2Cn,n−1Kn,n−1Mn,n2q2Cn,n−1Kn,nMn,n2q2Cn,nMn,n2000q2Cn,nKn,n−1Mn,n2q2Cn,nKn,nMn,n2000000000⋱00symmetryq2Kn,n−1Mn,n2q2Kn,n−1Kn,nMn,n2q1+q2Kn,nMn,n20,(20)R=r1r2⋱rm,where n is the degree of freedom of the structure. C, K, and M with subscripts represent the corresponding elements of the damping, stiffness, and mass matrices of the structure.


#### 3.1.2. Kalman Function Design

The index of control performance of the Kalman function is defined as(21)Je=limx⟶∞1T∫0Tzt−z^tTzt−z^tdt,where zt is the state variable vector defined in (2); z^t is the estimated state variable vector. The Kalman function is expressed as(22)z^˙t=Az^t+But+Keyt−y^ty^t=Cz^t,where A, B, C, and ut are defined in Equations (2)–(7); Ke is the gain matrix of the Kalman function. yt=x¨tTT is the measurement output, where x¨tTT is the acceleration vector. y^t is the estimated value of yt. Ke can be calculated as follows:(23)Ke=PeC0TRe−1.Using Riccati equation which is shown as(24)APe+PeAT−PeC0TRe−1C0Pe+Qe=0.In this study, we expressQe and Re as(25)Qe=αe,Re=βeI,where I is a unit matrix. αe and βe are coefficients of Qe and Re, respectively. Finally, we write the resultant controlling force of the LQG controller in equation (23) as(26)ftd=−GZ^t.


### 3.2. Hybrid Technique Combining WOA and LQG Controller

#### 3.2.1. Overview of WOA

The WOA is a intelligent optimization technique inspired by the hunting behavior of humpback whales. Each possible solution in the algorithm is denoted by codes carried by each whale. The solution corresponds to the whale’s position. Based on related literatures [26], the core mathematical model of WOA can be presented in two phases as follows. The pseudocode of this algorithm is described in Figure 1.Figure 1
Flowchart of WOA.(1) Searching and Encircling Prey. The mathematical model for searching prey is expressed as(27)Dw=Cw⋅Xrand−X,(28)Xt+1=Xrand−Aw⋅Dw,where t is the current iteration; Xrand is a random position vector selected from the current population. X is the position vector which needs updating in every iteration once a better solution is generated. Aw and Cw are coefficient respectively calculated according to(29)Aw=2a⋅r−a,Cw=2⋅r,where a has the range [2,0] with linear decreasement during the iteration. r is randomly chosen from the range [0,1].It is worth mentioning here that searching prey is conducted based on the above (27) and (28) on the condition that Aw≥1, which indicates that in this phase exploration is encouraged so that the algorithm can achieve a global search. If Aw≥1, shrinking mechanism for encircling prey will be performed. This is an exploitation phase. The current best solution can be exploited during this phase. Accordingly, update the position of each whale based on the equations as follows:(30)Dw′=Cw⋅X∗t−Xt,(31)Xt+1=X∗t−Aw⋅Dw′,where X∗ is the position vector of the best solution found so far.(2) Spiral Updating Position. To model the behavior of simultaneous encircling circularly and swimming spirally, it is assumed that the probabilities between these two movements to update the position of whales are both 50% [25–35]. This process can be mathematically modeled as(32)Xt+1=X∗t−Aw⋅Dw′ifp<0.5Dw″⋅ebl⋅cos2πl+X∗tifp<0.5,where p is randomly selected from the range 0∼1. Due to the probability 50%, the critical value of p for movements is 0.5. It should be noted that the latter movement is also an exploitation process. Dw″=X∗t−Xt is the distance from the ith whale to the best prey found so far. Constant b is used to define the spiral shape. l is a random value range from −1 to 1.


#### 3.2.2. Optimization Process

The WOA-based LQG controller is optimized according to the following steps.Step 1.
Determine theQe and Re for the gain matrix Ke of the Kalman function of the LQG controller.Step 2.
Define the function as [14](33)Obj=w×Obj1+1−w×Obj2,where(34)Obj1=maxt,ixitxuncmax,Obj2=maxt,ix¨aitx¨a,uncmax.
xi (t) is the relative displacement of the ith floor over the entire response. x¨ait is the absolute acceleration of ith floor over the entire response. Considering that the peak values of xi (t) and x¨ait may be negative, the absolute forms are adopted for both xi(t) and x¨ait. xuncmax and x¨a,uncmax represent the peak relative displacement and peak absolute acceleration, respectively, of the uncontrolled structure over the entire response. Obj1 and Obj2 are evaluation criteria and w is the weighting coefficient. In this study, for easy of calculation, the fitness function of WOA is equal to the objective function (see (32)). Besides, this is an optimization problem with constraints because it is necessary to limit the controlling force calculated by the WOA-optimized LQG within the range of the damping force (see Step 6).Step 3.
According to the characteristics of the controlled structure with MR dampers and the control objective, determine the number of inputs and outputs of the controller. Then deduce the structures of the weighting matricesQ and R of the optimal feedback gain.Step 4.
Determine the parameter boundaries and encoding strategy. In this study, we employ the real-value encoding strategy.Step 5.
Initialize the WOA: randomly generate initial locationX (0) = [X1, X2, ⋯, XN] for the population, and initialize the algorithm parameters, e.g. population size N and maximum iterative number Tmax.Step 6.
Calculate the optimal feedback gainG, feedback force fd, and the fitness value f(Xi) for each individual. Then find out the position with the best fitness value as the best position X∗. Let t = 1, go to Step 7. Notice that, according to Step 2, f(Xi) is taken as 1 if the controlling force is not within the damping force range; otherwise, f(Xi) is calculated based on equation (33).Step 7.
Begin iterative computation. Lett = t + 1, update a, A, C, l, and p;Step 8.
Ifp<0.5, the specific operations are as follows:
IfAw≥1, conduct searching prey. That is, randomly determine a position Xrand in the range of current population, and update the position Xi according to (28);
IfAw<1, conduct encircling prey. That is, update the position Xi based on X∗ according to (30).
Ifp≥0.5, conduct spiral position updating according to equation (31).Step 9.
After updating the position of every whale, judge whether the value of each parameter is out of the predefined range. If the value is larger than the upper limit (or smaller than the lower limit), replace it with the upper (or lower) limit.Step 10.
Calculate the fitness values for the whole population. If the best whale in the new population is better than the that in the original population, replace the latter with the former. Otherwise, keep the latter unchanged.Step 11.
Record the current best position and its corresponding fitness value. Ift < Tmax, go back to Step 7; otherwise, go to Step 12;Step 12.
Output the best position, which is the optimal parameter set for the weighting matricesQ and R.
The above-mentioned optimization flowchart can be described in Figure2.Figure 2
Optimization flowchart of the WLQG control.


## 3.1. Design of LQG Control Strategy

In this study, the LQG controller is designed to calculate the desired controlling force. It consists of two parts, i.e. regulator and observer filter which can calculate feedback force and estimate the full states separately.


### 3.1.1. Regulator Design

For the sake of security of the structure and the occupants’ comfort, the structural responses including maximum displacement and acceleration are required to be reduced. Besides, the desired controlling force should be limited in a reasonable range. Therefore, in this study, the index of control performance of the LQG is defined as(15)J=limx⟶∞1T∫0Tq1xmax2+q2x¨max2+∑i=1mrifdidt,where q1 and q2 are the weighting coefficients of the maximum relative displacement xmax and maximum absolute acceleration x¨max, respectively. For the sake of convenience, the relative displacement and absolute acceleration are called displacement and acceleration in short, respectively, in the subsequent discussion. The resultant feedback force fd which is a diagonal matrix consisting of fdi, also called the active force in the semi-active control, is used as one of the inputs for ANFIS inverse model to predict the command voltage (see Section 4.2).ri is the weighting coefficient of the feedback force fdi. M is the total number of the feedback forces, which is equal to the number of MR dampers. It should be noted that q1, q2, and ri are the parameters to be determined by WOA.The core of the control is to design a regulator so as to minimize the performance indexJ, which is rewritten as(16)J=limx⟶∞1T∫0TzTQz+fdTRfddt,where Q is the weighting matrix of state variable vector z defined in (2); R is the weighting matrix of the resultant feedback force fd defined in (15).According to the optimal control law, the resultant controlling force is written as(17)fdt=−12R−1BTPZ=−GZt,where G is the feedback gain, and P can be calculated by the following Riccati equation:(18)-PA-ATP+PBR−1BTP−Q=0,where A and B are the state equation matrices defined in (2) and (6), respectively, both of which are determined by the structural parameters.Finally, to calculate the optimal feedback gainG, the key step is to determine its weighting matrices Q and R. Based on (15) and (16), these weighting matrices are deduced as(19)Q=0000000000000000000⋱0000000q2Cn,n−1Mn,n2q2Cn,n−1Cn,nMn,n2000q2Cn,n−1Kn,n−1Mn,n2q2Cn,n−1Kn,nMn,n2q2Cn,nMn,n2000q2Cn,nKn,n−1Mn,n2q2Cn,nKn,nMn,n2000000000⋱00symmetryq2Kn,n−1Mn,n2q2Kn,n−1Kn,nMn,n2q1+q2Kn,nMn,n20,(20)R=r1r2⋱rm,where n is the degree of freedom of the structure. C, K, and M with subscripts represent the corresponding elements of the damping, stiffness, and mass matrices of the structure.


### 3.1.2. Kalman Function Design

The index of control performance of the Kalman function is defined as(21)Je=limx⟶∞1T∫0Tzt−z^tTzt−z^tdt,where zt is the state variable vector defined in (2); z^t is the estimated state variable vector. The Kalman function is expressed as(22)z^˙t=Az^t+But+Keyt−y^ty^t=Cz^t,where A, B, C, and ut are defined in Equations (2)–(7); Ke is the gain matrix of the Kalman function. yt=x¨tTT is the measurement output, where x¨tTT is the acceleration vector. y^t is the estimated value of yt. Ke can be calculated as follows:(23)Ke=PeC0TRe−1.Using Riccati equation which is shown as(24)APe+PeAT−PeC0TRe−1C0Pe+Qe=0.In this study, we expressQe and Re as(25)Qe=αe,Re=βeI,where I is a unit matrix. αe and βe are coefficients of Qe and Re, respectively. Finally, we write the resultant controlling force of the LQG controller in equation (23) as(26)ftd=−GZ^t.


## 3.1.1. Regulator Design

For the sake of security of the structure and the occupants’ comfort, the structural responses including maximum displacement and acceleration are required to be reduced. Besides, the desired controlling force should be limited in a reasonable range. Therefore, in this study, the index of control performance of the LQG is defined as(15)J=limx⟶∞1T∫0Tq1xmax2+q2x¨max2+∑i=1mrifdidt,where q1 and q2 are the weighting coefficients of the maximum relative displacement xmax and maximum absolute acceleration x¨max, respectively. For the sake of convenience, the relative displacement and absolute acceleration are called displacement and acceleration in short, respectively, in the subsequent discussion. The resultant feedback force fd which is a diagonal matrix consisting of fdi, also called the active force in the semi-active control, is used as one of the inputs for ANFIS inverse model to predict the command voltage (see Section 4.2).ri is the weighting coefficient of the feedback force fdi. M is the total number of the feedback forces, which is equal to the number of MR dampers. It should be noted that q1, q2, and ri are the parameters to be determined by WOA.The core of the control is to design a regulator so as to minimize the performance indexJ, which is rewritten as(16)J=limx⟶∞1T∫0TzTQz+fdTRfddt,where Q is the weighting matrix of state variable vector z defined in (2); R is the weighting matrix of the resultant feedback force fd defined in (15).According to the optimal control law, the resultant controlling force is written as(17)fdt=−12R−1BTPZ=−GZt,where G is the feedback gain, and P can be calculated by the following Riccati equation:(18)-PA-ATP+PBR−1BTP−Q=0,where A and B are the state equation matrices defined in (2) and (6), respectively, both of which are determined by the structural parameters.Finally, to calculate the optimal feedback gainG, the key step is to determine its weighting matrices Q and R. Based on (15) and (16), these weighting matrices are deduced as(19)Q=0000000000000000000⋱0000000q2Cn,n−1Mn,n2q2Cn,n−1Cn,nMn,n2000q2Cn,n−1Kn,n−1Mn,n2q2Cn,n−1Kn,nMn,n2q2Cn,nMn,n2000q2Cn,nKn,n−1Mn,n2q2Cn,nKn,nMn,n2000000000⋱00symmetryq2Kn,n−1Mn,n2q2Kn,n−1Kn,nMn,n2q1+q2Kn,nMn,n20,(20)R=r1r2⋱rm,where n is the degree of freedom of the structure. C, K, and M with subscripts represent the corresponding elements of the damping, stiffness, and mass matrices of the structure.


## 3.1.2. Kalman Function Design

The index of control performance of the Kalman function is defined as(21)Je=limx⟶∞1T∫0Tzt−z^tTzt−z^tdt,where zt is the state variable vector defined in (2); z^t is the estimated state variable vector. The Kalman function is expressed as(22)z^˙t=Az^t+But+Keyt−y^ty^t=Cz^t,where A, B, C, and ut are defined in Equations (2)–(7); Ke is the gain matrix of the Kalman function. yt=x¨tTT is the measurement output, where x¨tTT is the acceleration vector. y^t is the estimated value of yt. Ke can be calculated as follows:(23)Ke=PeC0TRe−1.Using Riccati equation which is shown as(24)APe+PeAT−PeC0TRe−1C0Pe+Qe=0.In this study, we expressQe and Re as(25)Qe=αe,Re=βeI,where I is a unit matrix. αe and βe are coefficients of Qe and Re, respectively. Finally, we write the resultant controlling force of the LQG controller in equation (23) as(26)ftd=−GZ^t.


## 3.2. Hybrid Technique Combining WOA and LQG Controller

### 3.2.1. Overview of WOA

The WOA is a intelligent optimization technique inspired by the hunting behavior of humpback whales. Each possible solution in the algorithm is denoted by codes carried by each whale. The solution corresponds to the whale’s position. Based on related literatures [26], the core mathematical model of WOA can be presented in two phases as follows. The pseudocode of this algorithm is described in Figure 1.Figure 1
Flowchart of WOA.(1) Searching and Encircling Prey. The mathematical model for searching prey is expressed as(27)Dw=Cw⋅Xrand−X,(28)Xt+1=Xrand−Aw⋅Dw,where t is the current iteration; Xrand is a random position vector selected from the current population. X is the position vector which needs updating in every iteration once a better solution is generated. Aw and Cw are coefficient respectively calculated according to(29)Aw=2a⋅r−a,Cw=2⋅r,where a has the range [2,0] with linear decreasement during the iteration. r is randomly chosen from the range [0,1].It is worth mentioning here that searching prey is conducted based on the above (27) and (28) on the condition that Aw≥1, which indicates that in this phase exploration is encouraged so that the algorithm can achieve a global search. If Aw≥1, shrinking mechanism for encircling prey will be performed. This is an exploitation phase. The current best solution can be exploited during this phase. Accordingly, update the position of each whale based on the equations as follows:(30)Dw′=Cw⋅X∗t−Xt,(31)Xt+1=X∗t−Aw⋅Dw′,where X∗ is the position vector of the best solution found so far.(2) Spiral Updating Position. To model the behavior of simultaneous encircling circularly and swimming spirally, it is assumed that the probabilities between these two movements to update the position of whales are both 50% [25–35]. This process can be mathematically modeled as(32)Xt+1=X∗t−Aw⋅Dw′ifp<0.5Dw″⋅ebl⋅cos2πl+X∗tifp<0.5,where p is randomly selected from the range 0∼1. Due to the probability 50%, the critical value of p for movements is 0.5. It should be noted that the latter movement is also an exploitation process. Dw″=X∗t−Xt is the distance from the ith whale to the best prey found so far. Constant b is used to define the spiral shape. l is a random value range from −1 to 1.


### 3.2.2. Optimization Process

The WOA-based LQG controller is optimized according to the following steps.Step 1.
Determine theQe and Re for the gain matrix Ke of the Kalman function of the LQG controller.Step 2.
Define the function as [14](33)Obj=w×Obj1+1−w×Obj2,where(34)Obj1=maxt,ixitxuncmax,Obj2=maxt,ix¨aitx¨a,uncmax.
xi (t) is the relative displacement of the ith floor over the entire response. x¨ait is the absolute acceleration of ith floor over the entire response. Considering that the peak values of xi (t) and x¨ait may be negative, the absolute forms are adopted for both xi(t) and x¨ait. xuncmax and x¨a,uncmax represent the peak relative displacement and peak absolute acceleration, respectively, of the uncontrolled structure over the entire response. Obj1 and Obj2 are evaluation criteria and w is the weighting coefficient. In this study, for easy of calculation, the fitness function of WOA is equal to the objective function (see (32)). Besides, this is an optimization problem with constraints because it is necessary to limit the controlling force calculated by the WOA-optimized LQG within the range of the damping force (see Step 6).Step 3.
According to the characteristics of the controlled structure with MR dampers and the control objective, determine the number of inputs and outputs of the controller. Then deduce the structures of the weighting matricesQ and R of the optimal feedback gain.Step 4.
Determine the parameter boundaries and encoding strategy. In this study, we employ the real-value encoding strategy.Step 5.
Initialize the WOA: randomly generate initial locationX (0) = [X1, X2, ⋯, XN] for the population, and initialize the algorithm parameters, e.g. population size N and maximum iterative number Tmax.Step 6.
Calculate the optimal feedback gainG, feedback force fd, and the fitness value f(Xi) for each individual. Then find out the position with the best fitness value as the best position X∗. Let t = 1, go to Step 7. Notice that, according to Step 2, f(Xi) is taken as 1 if the controlling force is not within the damping force range; otherwise, f(Xi) is calculated based on equation (33).Step 7.
Begin iterative computation. Lett = t + 1, update a, A, C, l, and p;Step 8.
Ifp<0.5, the specific operations are as follows:
IfAw≥1, conduct searching prey. That is, randomly determine a position Xrand in the range of current population, and update the position Xi according to (28);
IfAw<1, conduct encircling prey. That is, update the position Xi based on X∗ according to (30).
Ifp≥0.5, conduct spiral position updating according to equation (31).Step 9.
After updating the position of every whale, judge whether the value of each parameter is out of the predefined range. If the value is larger than the upper limit (or smaller than the lower limit), replace it with the upper (or lower) limit.Step 10.
Calculate the fitness values for the whole population. If the best whale in the new population is better than the that in the original population, replace the latter with the former. Otherwise, keep the latter unchanged.Step 11.
Record the current best position and its corresponding fitness value. Ift < Tmax, go back to Step 7; otherwise, go to Step 12;Step 12.
Output the best position, which is the optimal parameter set for the weighting matricesQ and R.
The above-mentioned optimization flowchart can be described in Figure2.Figure 2
Optimization flowchart of the WLQG control.


## 3.2.1. Overview of WOA

The WOA is a intelligent optimization technique inspired by the hunting behavior of humpback whales. Each possible solution in the algorithm is denoted by codes carried by each whale. The solution corresponds to the whale’s position. Based on related literatures [26], the core mathematical model of WOA can be presented in two phases as follows. The pseudocode of this algorithm is described in Figure 1.Figure 1
Flowchart of WOA.(1) Searching and Encircling Prey. The mathematical model for searching prey is expressed as(27)Dw=Cw⋅Xrand−X,(28)Xt+1=Xrand−Aw⋅Dw,where t is the current iteration; Xrand is a random position vector selected from the current population. X is the position vector which needs updating in every iteration once a better solution is generated. Aw and Cw are coefficient respectively calculated according to(29)Aw=2a⋅r−a,Cw=2⋅r,where a has the range [2,0] with linear decreasement during the iteration. r is randomly chosen from the range [0,1].It is worth mentioning here that searching prey is conducted based on the above (27) and (28) on the condition that Aw≥1, which indicates that in this phase exploration is encouraged so that the algorithm can achieve a global search. If Aw≥1, shrinking mechanism for encircling prey will be performed. This is an exploitation phase. The current best solution can be exploited during this phase. Accordingly, update the position of each whale based on the equations as follows:(30)Dw′=Cw⋅X∗t−Xt,(31)Xt+1=X∗t−Aw⋅Dw′,where X∗ is the position vector of the best solution found so far.(2) Spiral Updating Position. To model the behavior of simultaneous encircling circularly and swimming spirally, it is assumed that the probabilities between these two movements to update the position of whales are both 50% [25–35]. This process can be mathematically modeled as(32)Xt+1=X∗t−Aw⋅Dw′ifp<0.5Dw″⋅ebl⋅cos2πl+X∗tifp<0.5,where p is randomly selected from the range 0∼1. Due to the probability 50%, the critical value of p for movements is 0.5. It should be noted that the latter movement is also an exploitation process. Dw″=X∗t−Xt is the distance from the ith whale to the best prey found so far. Constant b is used to define the spiral shape. l is a random value range from −1 to 1.


## 3.2.2. Optimization Process

The WOA-based LQG controller is optimized according to the following steps.Step 1.
Determine theQe and Re for the gain matrix Ke of the Kalman function of the LQG controller.Step 2.
Define the function as [14](33)Obj=w×Obj1+1−w×Obj2,where(34)Obj1=maxt,ixitxuncmax,Obj2=maxt,ix¨aitx¨a,uncmax.
xi (t) is the relative displacement of the ith floor over the entire response. x¨ait is the absolute acceleration of ith floor over the entire response. Considering that the peak values of xi (t) and x¨ait may be negative, the absolute forms are adopted for both xi(t) and x¨ait. xuncmax and x¨a,uncmax represent the peak relative displacement and peak absolute acceleration, respectively, of the uncontrolled structure over the entire response. Obj1 and Obj2 are evaluation criteria and w is the weighting coefficient. In this study, for easy of calculation, the fitness function of WOA is equal to the objective function (see (32)). Besides, this is an optimization problem with constraints because it is necessary to limit the controlling force calculated by the WOA-optimized LQG within the range of the damping force (see Step 6).Step 3.
According to the characteristics of the controlled structure with MR dampers and the control objective, determine the number of inputs and outputs of the controller. Then deduce the structures of the weighting matricesQ and R of the optimal feedback gain.Step 4.
Determine the parameter boundaries and encoding strategy. In this study, we employ the real-value encoding strategy.Step 5.
Initialize the WOA: randomly generate initial locationX (0) = [X1, X2, ⋯, XN] for the population, and initialize the algorithm parameters, e.g. population size N and maximum iterative number Tmax.Step 6.
Calculate the optimal feedback gainG, feedback force fd, and the fitness value f(Xi) for each individual. Then find out the position with the best fitness value as the best position X∗. Let t = 1, go to Step 7. Notice that, according to Step 2, f(Xi) is taken as 1 if the controlling force is not within the damping force range; otherwise, f(Xi) is calculated based on equation (33).Step 7.
Begin iterative computation. Lett = t + 1, update a, A, C, l, and p;Step 8.
Ifp<0.5, the specific operations are as follows:
IfAw≥1, conduct searching prey. That is, randomly determine a position Xrand in the range of current population, and update the position Xi according to (28);
IfAw<1, conduct encircling prey. That is, update the position Xi based on X∗ according to (30).
Ifp≥0.5, conduct spiral position updating according to equation (31).Step 9.
After updating the position of every whale, judge whether the value of each parameter is out of the predefined range. If the value is larger than the upper limit (or smaller than the lower limit), replace it with the upper (or lower) limit.Step 10.
Calculate the fitness values for the whole population. If the best whale in the new population is better than the that in the original population, replace the latter with the former. Otherwise, keep the latter unchanged.Step 11.
Record the current best position and its corresponding fitness value. Ift < Tmax, go back to Step 7; otherwise, go to Step 12;Step 12.
Output the best position, which is the optimal parameter set for the weighting matricesQ and R.
The above-mentioned optimization flowchart can be described in Figure2.Figure 2
Optimization flowchart of the WLQG control.


## 4. Inverse Modeling the MR Damper with ANFIS Technique

### 4.1. ANFIS Architecture

ANFIS can accurately map the highly nonlinear relationship between inputs and outputs. Figure3 depicts an ANFIS with a two-input two-rule structure consisting of five layers [38].Figure 3
ANFIS with a two-input-two-rule structure [38].Layer 1.
The function of this layer is to fuzzify the input signals. Every nodei represents an output defined by(35)O1,i=μAix,i=1,2,O1,i=μBi−2y,i=3,4,where x and y are inputs to the nodes. Ai and Bi-2 denote linguistic labels. μA (or μB) is the membership function of a fuzzy set.Layer 2.
This layer is used to calculate the firing strength for each fuzzy rule, which is equal to the product of all the output signals of layer 1:(36)O2,i=wi=μAix×μBiy,i=1,2.Layer 3.
The function of this layer is to normalize the firing strength of each fuzzy rule by the following equation(37)O3,i=w¯i=wiw1+w2,i=1,2.Layer 4.
The function of this layer is to calculate the output of each rule, which can be expressed as(38)O4,i=w¯ifi=w¯ipix+qiy+ri,i=1,2,where {pi,qi,ri} is a set of adaptive rule parameters.Layer 5.
The function of this layer is to compute the overall output:(39)O5,i=∑i=12w¯ifi=∑i=12wifi∑i=12wi.


### 4.2. Inverse Modeling Strategy

Figure4 illustrates the MR damper inverse modeling using ANFIS technique. Given the displacement x(k), velocity x˙ (k), and target voltage u(k), the forward model built with a phenomenological model [36] is used to calculate the target damping force f(k). ANFIS with a four-input-one-output system is used to develop the MR damper inverse model. These inputs are x(k), x˙ (k), f(k), and the previous voltage u(k−1). Predicted voltage u^ (k) is the output of the system. According to the ANFIS training data, the inverse model will be optimally achieved through a hybrid training algorithm combining gradient descent and least square estimation methods [38]. The training process will terminate until the smallest root mean square error between the target and predicted voltages is attained or the iteration number is attained.Figure 4
The MR damper inverse modeling using ANFIS technique [39].Based on the predicted voltageu^ (k), the damping force f^ (k) can be computed through forward model for tracking the target damping force f(k). As shown in Figures 5–8 in Section 5, the effectiveness of the trained inverse model is validated in terms of comparison between f^ (k) and f(k). It should be noted that the inverse model is employed in a closed semiactive control system to calculate the MR damper input voltage. In this closed control system, the target damping force which is one of the ANFIS inputs is replaced by the feedback force fd of the LQG controller. In the meantime, the other three inputs of the ANFIS are replaced by the actual corresponding values in the closed control system.Figure 5
Predicted and target voltages (damping forces) for the training dataset.
(a)(b)Figure 6
Predicted and target voltages (damping forces) for validation dataset I.
(a)(b)Figure 7
Predicted and target voltages (damping forces) for validation dataset II.
(a)(b)Figure 8
Predicted and target voltages (damping forces) for validation dataset III.
(a)(b)


## 4.1. ANFIS Architecture

ANFIS can accurately map the highly nonlinear relationship between inputs and outputs. Figure3 depicts an ANFIS with a two-input two-rule structure consisting of five layers [38].Figure 3
ANFIS with a two-input-two-rule structure [38].Layer 1.
The function of this layer is to fuzzify the input signals. Every nodei represents an output defined by(35)O1,i=μAix,i=1,2,O1,i=μBi−2y,i=3,4,where x and y are inputs to the nodes. Ai and Bi-2 denote linguistic labels. μA (or μB) is the membership function of a fuzzy set.Layer 2.
This layer is used to calculate the firing strength for each fuzzy rule, which is equal to the product of all the output signals of layer 1:(36)O2,i=wi=μAix×μBiy,i=1,2.Layer 3.
The function of this layer is to normalize the firing strength of each fuzzy rule by the following equation(37)O3,i=w¯i=wiw1+w2,i=1,2.Layer 4.
The function of this layer is to calculate the output of each rule, which can be expressed as(38)O4,i=w¯ifi=w¯ipix+qiy+ri,i=1,2,where {pi,qi,ri} is a set of adaptive rule parameters.Layer 5.
The function of this layer is to compute the overall output:(39)O5,i=∑i=12w¯ifi=∑i=12wifi∑i=12wi.


## 4.2. Inverse Modeling Strategy

Figure4 illustrates the MR damper inverse modeling using ANFIS technique. Given the displacement x(k), velocity x˙ (k), and target voltage u(k), the forward model built with a phenomenological model [36] is used to calculate the target damping force f(k). ANFIS with a four-input-one-output system is used to develop the MR damper inverse model. These inputs are x(k), x˙ (k), f(k), and the previous voltage u(k−1). Predicted voltage u^ (k) is the output of the system. According to the ANFIS training data, the inverse model will be optimally achieved through a hybrid training algorithm combining gradient descent and least square estimation methods [38]. The training process will terminate until the smallest root mean square error between the target and predicted voltages is attained or the iteration number is attained.Figure 4
The MR damper inverse modeling using ANFIS technique [39].Based on the predicted voltageu^ (k), the damping force f^ (k) can be computed through forward model for tracking the target damping force f(k). As shown in Figures 5–8 in Section 5, the effectiveness of the trained inverse model is validated in terms of comparison between f^ (k) and f(k). It should be noted that the inverse model is employed in a closed semiactive control system to calculate the MR damper input voltage. In this closed control system, the target damping force which is one of the ANFIS inputs is replaced by the feedback force fd of the LQG controller. In the meantime, the other three inputs of the ANFIS are replaced by the actual corresponding values in the closed control system.Figure 5
Predicted and target voltages (damping forces) for the training dataset.
(a)(b)Figure 6
Predicted and target voltages (damping forces) for validation dataset I.
(a)(b)Figure 7
Predicted and target voltages (damping forces) for validation dataset II.
(a)(b)Figure 8
Predicted and target voltages (damping forces) for validation dataset III.
(a)(b)


## 5. Numerical Study

A 10-storey steel-frame structure installed with MR dampers with the maximum force 1000 kN and the maximum command voltage 9 V is considered for the numerical study [40]. The structural parameters are listed in Table 1. Considering that the 1940 El Centro earthquake is a far-field historical ground motion record usually chosen as the seismic excitation for the optimization design of the semi-active control strategy with MR dampers [17, 36, 37], in this study the building is subjected to the excitation of first 30 seconds of the 1940 El Centro earthquake with the maximum ground acceleration 3.417 m/s2. Considering that compared with the lower floors, the upper floors have larger uncontrolled relative displacement and absolute acceleration responses (see Figures 9 and 10), three MR dampers are installed respectively at the top three floors. To evaluate the performance of the proposed strategy, different control methods are used for comparison in this section.Table 1
Structural parameters of the building [40].
ParametersValuesStructural massmim1 = 2.15 × 105 kgmi = 2.01 × 105 kg, i = 2,3, …, 7mi = 2.03 × 105 kg, i = 8,9,10Structural stiffnesskik1 = 4.679 × 108 N/m, k2 = 4.76 × 108 N/mk3 = 4.68 × 108 N/mki = 4.5 × 108 N/m, i = 4,5,6,7ki = 4.37 × 108 N/m, i = 8,9,10Structural dampingcic1 = 1.676 × 106 Ns/m, c2 = 1.648 × 106 Ns/mc3 = 1.585 × 106 Ns/m, c4 = 1.585 × 106 Ns/mc5 = 1.539 × 106 Ns/m, c6 = 1.539 × 106 Ns/mc7 = 1.539 × 106 Ns/m, c8 = 1.539 × 106 Ns/mc9 = 1.099 × 106 Ns/m, c10 = 1.146 × 106 Ns/mFigure 9
Errors between predicted and target damping forces for validation datasets. (a) Error for validation dataset I. (b) Error for validation dataset II. (c) Error for validation dataset III.
(a)(b)(c)Figure 10
Comparison of the interstory drifts of all the floors among different control methods.


### 5.1. The Optimal WLQG Controller and Metaheuristic Algorithm Comparison

Since there are three MR dampers installed for the structure, the WOA-based LQG (here known as WLQG) controller is a 13-input-3-output system. Specifically, the inputs of the Kalman function of the controller formulated in (21) include three feedback forces and ten accelerations, and the outputs of the LQG controller include three feedback forces. The weighting matrices Q and R in equations (19) and (20) have dimensions of 10 × 10 and 3 × 3, respectively. In equation (33), w  = 0.1 is adopted as the weight value in the subsequent studies. In order to guarantee the good ability of mitigating the seismic responses, the value ranges of the parameters are chosen as below(40)q1,q2∈10−2,106r1,r2,r3∈10−6,103.Other three meta-heuristic algorithms, namely, GA, DE, and ABC are used for comparison. The controlling parameters of these three algorithms are listed in Table2. For a fair comparison, the population size and the maximum iterations of all these compared algorithms are set to be 30 and 120, respectively. Since these algorithms are stochastic methods, in this study every algorithm is executed for 20 independent runs having different initial population for obtaining average results [41, 42]. The convergence characteristics of the best solutions among 20 runs for these algorithms are shown in Figure 11. The figure shows that among these four algorithms WOA can achieve the smallest objective value and converge at the least number of iteration, indicating that WOA can obtain the best convergence accuracy and speed.Table 2
Initial parameters of the meta-heuristic algorithms.
AlgorithmsParameterValueGACrossover probability0.7Mutation probability0.02Selection mechanismStochastic universal samplingDEScaling factor for mutation0.4Crossover probability0.1ABCControl parameter to abandon the food source5Figure 11
Convergence curves of the best solutions obtained with different meta-heuristic algorithms.For a more specific comparison, the results in terms of the best, mean, and worst objective values and standard deviation of the objective values with different algorithms are presented in Table3. Boxplots in Figure 12 are utilized to visually determine the most effective algorithm. Table 3 shows that WOA outperforms its competitors in terms of all the indexes. Specifically, that WOA has the best convergence accuracy is further demonstrated due to its lowest mean value. Besides, this algorithm has the best robustness due to its lowest standard deviation. We can observe from Figure 12 that the smallest and median objective values obtained with WOA are both obviously smaller than those obtained with other three algorithms. These results can justify the best performance of WOA for optimizing the LQG controller. The reason why WOA outperforms other three algorithms is probably that it has good ability to switch between exploration and exploitation phases for local optima avoidance and high convergence speed.Table 3
Performance comparison among different algorithms.
Objective valueWOAGADEABCBest0.286700.295730.293830.28840Mean0.288260.300650.297850.29070Worst0.290840.306210.305540.29249Standard deviation0.001330.002960.002760.00140Figure 12
Boxplots of different metaheuristic algorithms for the objective value.To check whether the performance difference between WOA and any of the other three algorithms is statistically significant, ANOVA test is carried out with a 5% significance level. The null hypothesis assumes that there is no significant difference between compared algorithm; whereas, the alternative hypothesis considers a significant difference between them. Table4 shows that all the values of “Prob (probability)>F” for the statistical analyses are smaller than 5%, which contradicts the null hypothesis. Therefore, the ANOVA test demonstrates that there is a statistically significant difference between the WOA and GA/DE/ABC, and this conclusion does not just occur by coincidence.Table 4
Analysis results of ANOVA.
SourceSSMSFProb > FWOA/GA0.001540.00154291.842.00648E − 19WOA/DE0.000920.00092196.361.36046E − 16WOA/ABC0.000065.98463E − 0532.041.66596E − 06These algorithms are programmed in “Matlab R2016a 64-bit” and implemented on “Windows 7” environment on a computer having an 8.0 GB RAM-based Intel Core i7-4790 processor (3.6 GHz). The average execution time for WOA, GA, DE, and ABC is 7.3127E + 02 s, 9.6855 E + 02 s, 1.3649 E + 03 s, and 5.3811 E + 03 s, respectively. The reason why the runtime for WOA is the shortest is probably that this algorithm has low complexity and few parameters to be adjusted. ABC requires the longest runtime, although both of its convergence accuracy and robustness are better than GA and DE as shown in Table 3. However, how to design the controller is not a real-time problem. Therefore, we more concern about the quality of optimal solution than the execution time.


### 5.2. The MR Damper Inverse Model

#### 5.2.1. Training the Inverse Model

The training sample should be wide enough in terms of frequency and amplitude to ensure the validity of the inverse model. Therefore, as for the adopted MR damper, Gaussian white-noise signal having amplitude 1 m and frequency between 0 and 9 Hz is adopted to generate the displacement. The voltage is obtained with the same signals with amplitude between 0–10 V and frequency between 0–9 Hz, although the maximum input voltage is 9 V. Second-order backward difference method is adopted to calculate the displacement signals. Finally, based on these signals, forward model is employed to calculate the damping force. The training dataset has the duration of 10 seconds and is sampled at 1000 Hz. Therefore, 10000 sets of data are produced.The comparison results for the training dataset are illustrated in Figure5. Figure 5(a) shows that the predicted voltage matches very with the target one. Figure 5(b) compares the predicted damping force and the target one for the training dataset. This figure validates that there is almost no mismatch between these two damping forces, which means that the predicted damping force can accurately track the target one.


#### 5.2.2. Validation of the Inverse Model

For validating the ANFIS inverse model, we generate three groups of validation datasets different from the training dataset in terms of frequency and amplitude. These validation datasets are collected for 20 seconds and sampled at 250 Hz. The displacement and voltage signals used to produce the validation datasets are listed in Table5.Table 5
Displacement and voltage signals for validation datasets.
Validation datasetDisplacement (m)Voltage (V)IChip (0–3) Hz with amplitude 0.3Pulse with amplitude 5II0.5sin(4πt)7IIIChip (0–5) Hz with amplitude 0.86 + sin(3πt)Figures6–8 show the comparison results for the validation datasets. As shown in Figure 6(a) and Figure 7(a), it is obvious that the coincidence between these two voltages for the validation datasets is not as high as that for the training datasets. Errors between predicted and target damping forces for these three validation datasets are illustrated in Figures 11. Nevertheless, Figure 6(b)–8(b) indicate that the predicted damping force and the target one coincide quite well. In conclusion, the ANFIS inverse model possesses the capability to approximate the damper’s inverse dynamics.


#### 5.2.3. Control Results and Discussion

The proposed WLQG-ANFIS control strategy is compared with on-off control, LQG-ANFIS control, fuzzy control, and H∞-ANFIS control. In the LQG-ANFIS, Q and R are determined using(41)Q=aMK,R=bI,where K and M denote stiffness and mass matrices, respectively, as defined in (1). a and b are coefficients equal to 100 and 0.001, respectively. As for the fuzzy control, the accelerations of the top two floors are set to be inputs, and the outputs are the command voltages. Since this is a multi-input-multioutput fuzzy system, it is difficult to design the fuzzy inference rule based on expert experience; therefore in this study the fuzzy inference rule is tuned by GA [14]. In the H∞-ANFIS control, the top floor acceleration is used as the input of the H∞ controller. The inverse model adopted is the same as that in the WLQG-ANFIS control. The weighting function parameters of the H∞ controller are designed with GA by reference to [39].The peak structural responses including displacement, interstory drift and acceleration responses achieved by different control strategies are illustrated in Figures13, 10, and 14. The results demonstrate that these five control methods can decrease all the peak structural responses of all the floors except the peak accelerations of the two lowest floors even though there is not any MR damper installed on the lower seven floors. The proposed WLQG-ANFIS method shows superiority over other four control strategies in controlling all the peak responses even though minimizing the inter-story drifts is not the optimization objective for the WLQG controller design. Fuzzy control and on-off control outperform LQG-ANFIS and H∞-ANFIS in terms of all the upper floor peak responses. H∞-ANFIS performs slightly better than LQG-ANFIS for most of the peak responses.Figure 13
Comparison of the peak displacements of all the floors among different control methods.Figure 14
Comparison of peak accelerations of all the floors among different control methods.A set of criteria is used to evaluate the performance of these control methods which are defined as follows:(42)Δ1,i=xi,maxuc−xi,maxcxi,maxuc×100%,Δ2,i=di,maxuc−di,maxcdi,maxuc×100%,Δ3,i=x¨i,maxuc−x¨i,maxcx¨i,maxuc×100%.where xi,maxuc, di,maxuc, and x¨i,maxuc are the uncontrolled peak displacement, interstory drift and acceleration, respectively, of the ith floor. xi,maxc, di,maxc, and x¨i,maxc are the controlled peak displacement, interstory drift and acceleration respectively of the ith floor. These three kinds of performance indexes for different control methods are compared in Figures 15–17, respectively. Among these control methods, it is the WLQG-ANFIS that can achieve all the best quantities except Δ31. In other words, WLQG-ANFIS performs best in terms of all the peak displacement, interstory drift and acceleration responses of all the floors except the peak acceleration of the 1st floor. It can be concluded that compared with the other four control strategies, the proposed WLQG-ANFIS has more excellent comprehensive performance in reducing structural responses.Figure 15
Performance indexes for five control strategies for all the displacements.Figure 16
Performance indexes for five control strategies for all the interstory drifts.Figure 17
Performance indexes for five control strategies for all the accelerations.Table6 shows the comparison results of the maximum displacement xmax, interstory drift dmax, acceleration x¨max, and damping force Fmax among different control methods. In the table, reduction ratios of the responses with and without control are listed in the parentheses. All the highest reduction ratios in bold font are achieved by using the WLQG-ANFIS control strategy. Besides, a smaller maximum damping force is required by this control strategy compared with the other control methods except the fuzzy control, meaning that the generated damping force with the WLQG-ANFIS control can be fully utilized to reduce the seismic vibration. Fuzzy control and on-off control perform better than LQG-ANFIS and H∞-ANFIS in terms of all the maximum responses. Generally, the WLQG-ANFIS control can achieve better comprehensive performance than other four control strategies.Table 6
Maximum responses (reduction ratios) and damping force with different control strategies.
Max. responseUnc.On-offLQG-ANFISFuzzyH∞-ANFISWLQG-ANFISxmax (m)0.28490.1151 (60%)0.1362 (52%)0.1157 (59%)0.1306 (54%)0.1037(64%)dmax (m)0.04150.0216 (48%)0.0199 (52%)0.0165 (60%)0.0190 (54%)0.0147(65%)x¨max (m/s2)15.8197.8082 (51%)8.3342 (47%)7.7941 (51%)8.2258 (48%)7.5078 (53%)Fmax (kN)534460335428416Unc. Indicates without control.As the top floor has the largest uncontrolled displacement and acceleration responses among all the floors, Figures18 and 19 depict the displacement time history and acceleration time history of the top floor obtained with the proposed control, fuzzy control, and on-off control. The results clearly show that the WLQG-ANFIS method is superior to other two methods in both responses.Figure 18
The top floor displacement responses with different control methods.Figure 19
The top floor acceleration responses with different control methods.Next, to investigate the robustness of the WLQG-ANFIS under the untrained earthquake ground excitations, the first 30 seconds of El Centro, Kobe, and Wenchuan ground acceleration records with the scaled maximum values of 8.54 m/s2, 5.72 m/s2, and 9.15 m/s2, respectively, are used as the excitations. Figures 20–22 illustrate these three scaled earthquake accelerations. Figures 23–25 illustrate the peak responses of all the floors with and without control under the scaled El Centro, Kobe, and Wenchuan earthquakes, respectively. The comparison of the maximum responses with the proposed control and without control subjected to these earthquake excitations are shown in Table 7. These results verify that the WLQG-ANFIS control can still perform satisfactorily even though the structure undergoes untrained earthquake excitations.Figure 20
The scaled El-centro earthquake accelerations.Figure 21
The scaled Kobe earthquake acceleration.Figure 22
The scaled Wenchuan earthquake acceleration.Figure 23
Peak responses due to scaled El-centro earthquake.
(a)(b)(c)Figure 24
Peak responses due to scaled Kobe earthquake.
(a)(b)(c)Figure 25
Peak responses due to scaled Wenchuan earthquake.
(a)(b)(c)Table 7
Maximum responses (reduction ratios) under varied seismic excitations with and without control.
Earthquake excitationControl strategyMaximum responsesxmax (m)dmax (m)x¨max (m/s2)Scaled el-centroUnc.0.71230.103839.5471Ctr.0.2711 (62%)0.0384 (63%)18.6426 (53%)Scaled kobeUnc.0.47590.069027.1337Ctr.0.2940 (38%)0.0409 (41%)16.6455 (39%)Scaled wenchuanUnc.0.20660.030120.0012Ctr.0.1007 (51%)0.0133 (56%)12.6535 (37%)Unc indicates without control; Ctr indicates the proposed control.In addition, to investigate the robustness of the WLQG-ANFIS method with regard to uncertainties in structural parameters owing to damage, structures having +15% and -15% stiffness variations for each floor are respectively considered. The peak responses of all the floors with and without control considering stiffness variations are illustrated in Figures26 and 27, respectively. The comparison of the maximum responses with the proposed control and without control for these two stiffness variation cases is shown in Table 8. These results demonstrate that the WLQG-ANFIS method can reduce the structural responses evidently although structural damages exist.Figure 26
Peak responses with consideration of 15% stiffness variation.
(a)(b)(c)Figure 27
Peak responses with consideration of −15% stiffness variation.
(a)(b)(c)Table 8
Maximum responses (reduction ratios) in different stiffness variation cases with the proposed control.
Stiffness variationControl strategyMaximum responsesxmax (m)dmax (m)x¨max (m/s2)+15%Unc.0.21010.033212.3863Ctr.0.1006 (52%)0.0131 (60%)8.8552 (29%)−15%Unc.0.26000.036914.6520Ctr.0.1169 (55%)0.0178 (52%)7.5895 (48%)Unc indicates without control; Ctr indicates the proposed control.


## 5.1. The Optimal WLQG Controller and Metaheuristic Algorithm Comparison

Since there are three MR dampers installed for the structure, the WOA-based LQG (here known as WLQG) controller is a 13-input-3-output system. Specifically, the inputs of the Kalman function of the controller formulated in (21) include three feedback forces and ten accelerations, and the outputs of the LQG controller include three feedback forces. The weighting matrices Q and R in equations (19) and (20) have dimensions of 10 × 10 and 3 × 3, respectively. In equation (33), w  = 0.1 is adopted as the weight value in the subsequent studies. In order to guarantee the good ability of mitigating the seismic responses, the value ranges of the parameters are chosen as below(40)q1,q2∈10−2,106r1,r2,r3∈10−6,103.Other three meta-heuristic algorithms, namely, GA, DE, and ABC are used for comparison. The controlling parameters of these three algorithms are listed in Table2. For a fair comparison, the population size and the maximum iterations of all these compared algorithms are set to be 30 and 120, respectively. Since these algorithms are stochastic methods, in this study every algorithm is executed for 20 independent runs having different initial population for obtaining average results [41, 42]. The convergence characteristics of the best solutions among 20 runs for these algorithms are shown in Figure 11. The figure shows that among these four algorithms WOA can achieve the smallest objective value and converge at the least number of iteration, indicating that WOA can obtain the best convergence accuracy and speed.Table 2
Initial parameters of the meta-heuristic algorithms.
AlgorithmsParameterValueGACrossover probability0.7Mutation probability0.02Selection mechanismStochastic universal samplingDEScaling factor for mutation0.4Crossover probability0.1ABCControl parameter to abandon the food source5Figure 11
Convergence curves of the best solutions obtained with different meta-heuristic algorithms.For a more specific comparison, the results in terms of the best, mean, and worst objective values and standard deviation of the objective values with different algorithms are presented in Table3. Boxplots in Figure 12 are utilized to visually determine the most effective algorithm. Table 3 shows that WOA outperforms its competitors in terms of all the indexes. Specifically, that WOA has the best convergence accuracy is further demonstrated due to its lowest mean value. Besides, this algorithm has the best robustness due to its lowest standard deviation. We can observe from Figure 12 that the smallest and median objective values obtained with WOA are both obviously smaller than those obtained with other three algorithms. These results can justify the best performance of WOA for optimizing the LQG controller. The reason why WOA outperforms other three algorithms is probably that it has good ability to switch between exploration and exploitation phases for local optima avoidance and high convergence speed.Table 3
Performance comparison among different algorithms.
Objective valueWOAGADEABCBest0.286700.295730.293830.28840Mean0.288260.300650.297850.29070Worst0.290840.306210.305540.29249Standard deviation0.001330.002960.002760.00140Figure 12
Boxplots of different metaheuristic algorithms for the objective value.To check whether the performance difference between WOA and any of the other three algorithms is statistically significant, ANOVA test is carried out with a 5% significance level. The null hypothesis assumes that there is no significant difference between compared algorithm; whereas, the alternative hypothesis considers a significant difference between them. Table4 shows that all the values of “Prob (probability)>F” for the statistical analyses are smaller than 5%, which contradicts the null hypothesis. Therefore, the ANOVA test demonstrates that there is a statistically significant difference between the WOA and GA/DE/ABC, and this conclusion does not just occur by coincidence.Table 4
Analysis results of ANOVA.
SourceSSMSFProb > FWOA/GA0.001540.00154291.842.00648E − 19WOA/DE0.000920.00092196.361.36046E − 16WOA/ABC0.000065.98463E − 0532.041.66596E − 06These algorithms are programmed in “Matlab R2016a 64-bit” and implemented on “Windows 7” environment on a computer having an 8.0 GB RAM-based Intel Core i7-4790 processor (3.6 GHz). The average execution time for WOA, GA, DE, and ABC is 7.3127E + 02 s, 9.6855 E + 02 s, 1.3649 E + 03 s, and 5.3811 E + 03 s, respectively. The reason why the runtime for WOA is the shortest is probably that this algorithm has low complexity and few parameters to be adjusted. ABC requires the longest runtime, although both of its convergence accuracy and robustness are better than GA and DE as shown in Table 3. However, how to design the controller is not a real-time problem. Therefore, we more concern about the quality of optimal solution than the execution time.


## 5.2. The MR Damper Inverse Model

### 5.2.1. Training the Inverse Model

The training sample should be wide enough in terms of frequency and amplitude to ensure the validity of the inverse model. Therefore, as for the adopted MR damper, Gaussian white-noise signal having amplitude 1 m and frequency between 0 and 9 Hz is adopted to generate the displacement. The voltage is obtained with the same signals with amplitude between 0–10 V and frequency between 0–9 Hz, although the maximum input voltage is 9 V. Second-order backward difference method is adopted to calculate the displacement signals. Finally, based on these signals, forward model is employed to calculate the damping force. The training dataset has the duration of 10 seconds and is sampled at 1000 Hz. Therefore, 10000 sets of data are produced.The comparison results for the training dataset are illustrated in Figure5. Figure 5(a) shows that the predicted voltage matches very with the target one. Figure 5(b) compares the predicted damping force and the target one for the training dataset. This figure validates that there is almost no mismatch between these two damping forces, which means that the predicted damping force can accurately track the target one.


### 5.2.2. Validation of the Inverse Model

For validating the ANFIS inverse model, we generate three groups of validation datasets different from the training dataset in terms of frequency and amplitude. These validation datasets are collected for 20 seconds and sampled at 250 Hz. The displacement and voltage signals used to produce the validation datasets are listed in Table5.Table 5
Displacement and voltage signals for validation datasets.
Validation datasetDisplacement (m)Voltage (V)IChip (0–3) Hz with amplitude 0.3Pulse with amplitude 5II0.5sin(4πt)7IIIChip (0–5) Hz with amplitude 0.86 + sin(3πt)Figures6–8 show the comparison results for the validation datasets. As shown in Figure 6(a) and Figure 7(a), it is obvious that the coincidence between these two voltages for the validation datasets is not as high as that for the training datasets. Errors between predicted and target damping forces for these three validation datasets are illustrated in Figures 11. Nevertheless, Figure 6(b)–8(b) indicate that the predicted damping force and the target one coincide quite well. In conclusion, the ANFIS inverse model possesses the capability to approximate the damper’s inverse dynamics.


### 5.2.3. Control Results and Discussion

The proposed WLQG-ANFIS control strategy is compared with on-off control, LQG-ANFIS control, fuzzy control, and H∞-ANFIS control. In the LQG-ANFIS, Q and R are determined using(41)Q=aMK,R=bI,where K and M denote stiffness and mass matrices, respectively, as defined in (1). a and b are coefficients equal to 100 and 0.001, respectively. As for the fuzzy control, the accelerations of the top two floors are set to be inputs, and the outputs are the command voltages. Since this is a multi-input-multioutput fuzzy system, it is difficult to design the fuzzy inference rule based on expert experience; therefore in this study the fuzzy inference rule is tuned by GA [14]. In the H∞-ANFIS control, the top floor acceleration is used as the input of the H∞ controller. The inverse model adopted is the same as that in the WLQG-ANFIS control. The weighting function parameters of the H∞ controller are designed with GA by reference to [39].The peak structural responses including displacement, interstory drift and acceleration responses achieved by different control strategies are illustrated in Figures13, 10, and 14. The results demonstrate that these five control methods can decrease all the peak structural responses of all the floors except the peak accelerations of the two lowest floors even though there is not any MR damper installed on the lower seven floors. The proposed WLQG-ANFIS method shows superiority over other four control strategies in controlling all the peak responses even though minimizing the inter-story drifts is not the optimization objective for the WLQG controller design. Fuzzy control and on-off control outperform LQG-ANFIS and H∞-ANFIS in terms of all the upper floor peak responses. H∞-ANFIS performs slightly better than LQG-ANFIS for most of the peak responses.Figure 13
Comparison of the peak displacements of all the floors among different control methods.Figure 14
Comparison of peak accelerations of all the floors among different control methods.A set of criteria is used to evaluate the performance of these control methods which are defined as follows:(42)Δ1,i=xi,maxuc−xi,maxcxi,maxuc×100%,Δ2,i=di,maxuc−di,maxcdi,maxuc×100%,Δ3,i=x¨i,maxuc−x¨i,maxcx¨i,maxuc×100%.where xi,maxuc, di,maxuc, and x¨i,maxuc are the uncontrolled peak displacement, interstory drift and acceleration, respectively, of the ith floor. xi,maxc, di,maxc, and x¨i,maxc are the controlled peak displacement, interstory drift and acceleration respectively of the ith floor. These three kinds of performance indexes for different control methods are compared in Figures 15–17, respectively. Among these control methods, it is the WLQG-ANFIS that can achieve all the best quantities except Δ31. In other words, WLQG-ANFIS performs best in terms of all the peak displacement, interstory drift and acceleration responses of all the floors except the peak acceleration of the 1st floor. It can be concluded that compared with the other four control strategies, the proposed WLQG-ANFIS has more excellent comprehensive performance in reducing structural responses.Figure 15
Performance indexes for five control strategies for all the displacements.Figure 16
Performance indexes for five control strategies for all the interstory drifts.Figure 17
Performance indexes for five control strategies for all the accelerations.Table6 shows the comparison results of the maximum displacement xmax, interstory drift dmax, acceleration x¨max, and damping force Fmax among different control methods. In the table, reduction ratios of the responses with and without control are listed in the parentheses. All the highest reduction ratios in bold font are achieved by using the WLQG-ANFIS control strategy. Besides, a smaller maximum damping force is required by this control strategy compared with the other control methods except the fuzzy control, meaning that the generated damping force with the WLQG-ANFIS control can be fully utilized to reduce the seismic vibration. Fuzzy control and on-off control perform better than LQG-ANFIS and H∞-ANFIS in terms of all the maximum responses. Generally, the WLQG-ANFIS control can achieve better comprehensive performance than other four control strategies.Table 6
Maximum responses (reduction ratios) and damping force with different control strategies.
Max. responseUnc.On-offLQG-ANFISFuzzyH∞-ANFISWLQG-ANFISxmax (m)0.28490.1151 (60%)0.1362 (52%)0.1157 (59%)0.1306 (54%)0.1037(64%)dmax (m)0.04150.0216 (48%)0.0199 (52%)0.0165 (60%)0.0190 (54%)0.0147(65%)x¨max (m/s2)15.8197.8082 (51%)8.3342 (47%)7.7941 (51%)8.2258 (48%)7.5078 (53%)Fmax (kN)534460335428416Unc. Indicates without control.As the top floor has the largest uncontrolled displacement and acceleration responses among all the floors, Figures18 and 19 depict the displacement time history and acceleration time history of the top floor obtained with the proposed control, fuzzy control, and on-off control. The results clearly show that the WLQG-ANFIS method is superior to other two methods in both responses.Figure 18
The top floor displacement responses with different control methods.Figure 19
The top floor acceleration responses with different control methods.Next, to investigate the robustness of the WLQG-ANFIS under the untrained earthquake ground excitations, the first 30 seconds of El Centro, Kobe, and Wenchuan ground acceleration records with the scaled maximum values of 8.54 m/s2, 5.72 m/s2, and 9.15 m/s2, respectively, are used as the excitations. Figures 20–22 illustrate these three scaled earthquake accelerations. Figures 23–25 illustrate the peak responses of all the floors with and without control under the scaled El Centro, Kobe, and Wenchuan earthquakes, respectively. The comparison of the maximum responses with the proposed control and without control subjected to these earthquake excitations are shown in Table 7. These results verify that the WLQG-ANFIS control can still perform satisfactorily even though the structure undergoes untrained earthquake excitations.Figure 20
The scaled El-centro earthquake accelerations.Figure 21
The scaled Kobe earthquake acceleration.Figure 22
The scaled Wenchuan earthquake acceleration.Figure 23
Peak responses due to scaled El-centro earthquake.
(a)(b)(c)Figure 24
Peak responses due to scaled Kobe earthquake.
(a)(b)(c)Figure 25
Peak responses due to scaled Wenchuan earthquake.
(a)(b)(c)Table 7
Maximum responses (reduction ratios) under varied seismic excitations with and without control.
Earthquake excitationControl strategyMaximum responsesxmax (m)dmax (m)x¨max (m/s2)Scaled el-centroUnc.0.71230.103839.5471Ctr.0.2711 (62%)0.0384 (63%)18.6426 (53%)Scaled kobeUnc.0.47590.069027.1337Ctr.0.2940 (38%)0.0409 (41%)16.6455 (39%)Scaled wenchuanUnc.0.20660.030120.0012Ctr.0.1007 (51%)0.0133 (56%)12.6535 (37%)Unc indicates without control; Ctr indicates the proposed control.In addition, to investigate the robustness of the WLQG-ANFIS method with regard to uncertainties in structural parameters owing to damage, structures having +15% and -15% stiffness variations for each floor are respectively considered. The peak responses of all the floors with and without control considering stiffness variations are illustrated in Figures26 and 27, respectively. The comparison of the maximum responses with the proposed control and without control for these two stiffness variation cases is shown in Table 8. These results demonstrate that the WLQG-ANFIS method can reduce the structural responses evidently although structural damages exist.Figure 26
Peak responses with consideration of 15% stiffness variation.
(a)(b)(c)Figure 27
Peak responses with consideration of −15% stiffness variation.
(a)(b)(c)Table 8
Maximum responses (reduction ratios) in different stiffness variation cases with the proposed control.
Stiffness variationControl strategyMaximum responsesxmax (m)dmax (m)x¨max (m/s2)+15%Unc.0.21010.033212.3863Ctr.0.1006 (52%)0.0131 (60%)8.8552 (29%)−15%Unc.0.26000.036914.6520Ctr.0.1169 (55%)0.0178 (52%)7.5895 (48%)Unc indicates without control; Ctr indicates the proposed control.


## 5.2.1. Training the Inverse Model

The training sample should be wide enough in terms of frequency and amplitude to ensure the validity of the inverse model. Therefore, as for the adopted MR damper, Gaussian white-noise signal having amplitude 1 m and frequency between 0 and 9 Hz is adopted to generate the displacement. The voltage is obtained with the same signals with amplitude between 0–10 V and frequency between 0–9 Hz, although the maximum input voltage is 9 V. Second-order backward difference method is adopted to calculate the displacement signals. Finally, based on these signals, forward model is employed to calculate the damping force. The training dataset has the duration of 10 seconds and is sampled at 1000 Hz. Therefore, 10000 sets of data are produced.The comparison results for the training dataset are illustrated in Figure5. Figure 5(a) shows that the predicted voltage matches very with the target one. Figure 5(b) compares the predicted damping force and the target one for the training dataset. This figure validates that there is almost no mismatch between these two damping forces, which means that the predicted damping force can accurately track the target one.


## 5.2.2. Validation of the Inverse Model

For validating the ANFIS inverse model, we generate three groups of validation datasets different from the training dataset in terms of frequency and amplitude. These validation datasets are collected for 20 seconds and sampled at 250 Hz. The displacement and voltage signals used to produce the validation datasets are listed in Table5.Table 5
Displacement and voltage signals for validation datasets.
Validation datasetDisplacement (m)Voltage (V)IChip (0–3) Hz with amplitude 0.3Pulse with amplitude 5II0.5sin(4πt)7IIIChip (0–5) Hz with amplitude 0.86 + sin(3πt)Figures6–8 show the comparison results for the validation datasets. As shown in Figure 6(a) and Figure 7(a), it is obvious that the coincidence between these two voltages for the validation datasets is not as high as that for the training datasets. Errors between predicted and target damping forces for these three validation datasets are illustrated in Figures 11. Nevertheless, Figure 6(b)–8(b) indicate that the predicted damping force and the target one coincide quite well. In conclusion, the ANFIS inverse model possesses the capability to approximate the damper’s inverse dynamics.


## 5.2.3. Control Results and Discussion

The proposed WLQG-ANFIS control strategy is compared with on-off control, LQG-ANFIS control, fuzzy control, and H∞-ANFIS control. In the LQG-ANFIS, Q and R are determined using(41)Q=aMK,R=bI,where K and M denote stiffness and mass matrices, respectively, as defined in (1). a and b are coefficients equal to 100 and 0.001, respectively. As for the fuzzy control, the accelerations of the top two floors are set to be inputs, and the outputs are the command voltages. Since this is a multi-input-multioutput fuzzy system, it is difficult to design the fuzzy inference rule based on expert experience; therefore in this study the fuzzy inference rule is tuned by GA [14]. In the H∞-ANFIS control, the top floor acceleration is used as the input of the H∞ controller. The inverse model adopted is the same as that in the WLQG-ANFIS control. The weighting function parameters of the H∞ controller are designed with GA by reference to [39].The peak structural responses including displacement, interstory drift and acceleration responses achieved by different control strategies are illustrated in Figures13, 10, and 14. The results demonstrate that these five control methods can decrease all the peak structural responses of all the floors except the peak accelerations of the two lowest floors even though there is not any MR damper installed on the lower seven floors. The proposed WLQG-ANFIS method shows superiority over other four control strategies in controlling all the peak responses even though minimizing the inter-story drifts is not the optimization objective for the WLQG controller design. Fuzzy control and on-off control outperform LQG-ANFIS and H∞-ANFIS in terms of all the upper floor peak responses. H∞-ANFIS performs slightly better than LQG-ANFIS for most of the peak responses.Figure 13
Comparison of the peak displacements of all the floors among different control methods.Figure 14
Comparison of peak accelerations of all the floors among different control methods.A set of criteria is used to evaluate the performance of these control methods which are defined as follows:(42)Δ1,i=xi,maxuc−xi,maxcxi,maxuc×100%,Δ2,i=di,maxuc−di,maxcdi,maxuc×100%,Δ3,i=x¨i,maxuc−x¨i,maxcx¨i,maxuc×100%.where xi,maxuc, di,maxuc, and x¨i,maxuc are the uncontrolled peak displacement, interstory drift and acceleration, respectively, of the ith floor. xi,maxc, di,maxc, and x¨i,maxc are the controlled peak displacement, interstory drift and acceleration respectively of the ith floor. These three kinds of performance indexes for different control methods are compared in Figures 15–17, respectively. Among these control methods, it is the WLQG-ANFIS that can achieve all the best quantities except Δ31. In other words, WLQG-ANFIS performs best in terms of all the peak displacement, interstory drift and acceleration responses of all the floors except the peak acceleration of the 1st floor. It can be concluded that compared with the other four control strategies, the proposed WLQG-ANFIS has more excellent comprehensive performance in reducing structural responses.Figure 15
Performance indexes for five control strategies for all the displacements.Figure 16
Performance indexes for five control strategies for all the interstory drifts.Figure 17
Performance indexes for five control strategies for all the accelerations.Table6 shows the comparison results of the maximum displacement xmax, interstory drift dmax, acceleration x¨max, and damping force Fmax among different control methods. In the table, reduction ratios of the responses with and without control are listed in the parentheses. All the highest reduction ratios in bold font are achieved by using the WLQG-ANFIS control strategy. Besides, a smaller maximum damping force is required by this control strategy compared with the other control methods except the fuzzy control, meaning that the generated damping force with the WLQG-ANFIS control can be fully utilized to reduce the seismic vibration. Fuzzy control and on-off control perform better than LQG-ANFIS and H∞-ANFIS in terms of all the maximum responses. Generally, the WLQG-ANFIS control can achieve better comprehensive performance than other four control strategies.Table 6
Maximum responses (reduction ratios) and damping force with different control strategies.
Max. responseUnc.On-offLQG-ANFISFuzzyH∞-ANFISWLQG-ANFISxmax (m)0.28490.1151 (60%)0.1362 (52%)0.1157 (59%)0.1306 (54%)0.1037(64%)dmax (m)0.04150.0216 (48%)0.0199 (52%)0.0165 (60%)0.0190 (54%)0.0147(65%)x¨max (m/s2)15.8197.8082 (51%)8.3342 (47%)7.7941 (51%)8.2258 (48%)7.5078 (53%)Fmax (kN)534460335428416Unc. Indicates without control.As the top floor has the largest uncontrolled displacement and acceleration responses among all the floors, Figures18 and 19 depict the displacement time history and acceleration time history of the top floor obtained with the proposed control, fuzzy control, and on-off control. The results clearly show that the WLQG-ANFIS method is superior to other two methods in both responses.Figure 18
The top floor displacement responses with different control methods.Figure 19
The top floor acceleration responses with different control methods.Next, to investigate the robustness of the WLQG-ANFIS under the untrained earthquake ground excitations, the first 30 seconds of El Centro, Kobe, and Wenchuan ground acceleration records with the scaled maximum values of 8.54 m/s2, 5.72 m/s2, and 9.15 m/s2, respectively, are used as the excitations. Figures 20–22 illustrate these three scaled earthquake accelerations. Figures 23–25 illustrate the peak responses of all the floors with and without control under the scaled El Centro, Kobe, and Wenchuan earthquakes, respectively. The comparison of the maximum responses with the proposed control and without control subjected to these earthquake excitations are shown in Table 7. These results verify that the WLQG-ANFIS control can still perform satisfactorily even though the structure undergoes untrained earthquake excitations.Figure 20
The scaled El-centro earthquake accelerations.Figure 21
The scaled Kobe earthquake acceleration.Figure 22
The scaled Wenchuan earthquake acceleration.Figure 23
Peak responses due to scaled El-centro earthquake.
(a)(b)(c)Figure 24
Peak responses due to scaled Kobe earthquake.
(a)(b)(c)Figure 25
Peak responses due to scaled Wenchuan earthquake.
(a)(b)(c)Table 7
Maximum responses (reduction ratios) under varied seismic excitations with and without control.
Earthquake excitationControl strategyMaximum responsesxmax (m)dmax (m)x¨max (m/s2)Scaled el-centroUnc.0.71230.103839.5471Ctr.0.2711 (62%)0.0384 (63%)18.6426 (53%)Scaled kobeUnc.0.47590.069027.1337Ctr.0.2940 (38%)0.0409 (41%)16.6455 (39%)Scaled wenchuanUnc.0.20660.030120.0012Ctr.0.1007 (51%)0.0133 (56%)12.6535 (37%)Unc indicates without control; Ctr indicates the proposed control.In addition, to investigate the robustness of the WLQG-ANFIS method with regard to uncertainties in structural parameters owing to damage, structures having +15% and -15% stiffness variations for each floor are respectively considered. The peak responses of all the floors with and without control considering stiffness variations are illustrated in Figures26 and 27, respectively. The comparison of the maximum responses with the proposed control and without control for these two stiffness variation cases is shown in Table 8. These results demonstrate that the WLQG-ANFIS method can reduce the structural responses evidently although structural damages exist.Figure 26
Peak responses with consideration of 15% stiffness variation.
(a)(b)(c)Figure 27
Peak responses with consideration of −15% stiffness variation.
(a)(b)(c)Table 8
Maximum responses (reduction ratios) in different stiffness variation cases with the proposed control.
Stiffness variationControl strategyMaximum responsesxmax (m)dmax (m)x¨max (m/s2)+15%Unc.0.21010.033212.3863Ctr.0.1006 (52%)0.0131 (60%)8.8552 (29%)−15%Unc.0.26000.036914.6520Ctr.0.1169 (55%)0.0178 (52%)7.5895 (48%)Unc indicates without control; Ctr indicates the proposed control.


## 6. Conclusion

An intelligent semiactive control method is presented in this study to mitigate the seismic vibration using MR dampers. In this semiactive control method, a parameter identification approach to optimize the LQG control using the WOA is proposed for the first time. The design of the LQG control is transformed into a constraint-based optimization problem related to the index of control performance. The WOA-based LQG controller and the MR damper ANFIS inverse modeling technique are combined for the semiactive control.The analysis results for a seismically-excited 10-storey structure installed with MR dampers are mainly concluded as follows: ANOVA test results statistically show that there is significant difference between the WOA and other three algorithms. The performance of the WOA surpasses that of GA, DE, and ABC in terms of convergence performance and robustness, and this conclusion is not just a matter of chance. The inverse dynamics of the MR damper can be characterized by the inverse model designed with the ANFIS technology, through which the predicted force can track the desired one very well. The proposed semiactive control method can evidently reduce the maximum displacement response and acceleration response with much low force requirement. Meanwhile, the interstory drifts of each floor can be effectively decreased. This control strategy exhibits more excellent comprehensive performance than on-off control, LQG-ANFIS control, fuzzy control, and H∞-ANFIS control. Besides, the proposed strategy can still satisfactorily mitigate structural vibration even considering the changes in structural parameters and seismic excitations.Efforts are currently underway to improve the WOA optimization of LQG control in order to attain better control results, and compare the algorithm with other more meta-heuristic algorithms. Experimental investigation to verify the proposed control scheme and other novel control schemes is also involved in the future research work. Besides, our future work also involves experimental investigation on the effect of varied temperature of MR fluid on the performance of MR damper and control results in practical application.



---

*Source: 4060660-2022-03-28.xml*