**1. Introduction**

Structural control systems have turned into a standard technology to improve the dynamic response of civil engineering structures subjected to dynamic actions, such as wind forces or earthquake loads [1]. These control systems can be classified into four major groups: passive, active, hybrid, and semiactive controllers. Passive systems are widely accepted by the engineering community because of their mechanical simplicity, low power requirements, and controllable force capacity [2]. Among passive systems, one of the most commonly used and tested devices throughout the years has been the tuned mass damper (TMD). This system consists of attaching an additional mass linked to the main structure using a spring and a viscous damper, which is optimally tuned to one of the fundamental vibration frequencies of the system in order to transfer energy among the vibrating modes by making the structure more flexible [3]. The original TMD formulation was first proposed by Frahm in [4] as a vibration absorber with no damping to control periodic resonance vibrations. Subsequently, Ormondroyd and Den Hartog [5] developed the theory of vibration absorbers, including viscous damping to the system to be effective under different frequencies of random vibrations. Conventional tuning methodologies

**Citation:** Lara-Valencia, L.A.; Caicedo, D.; Valencia-Gonzalez, Y. A Novel Whale Optimization Algorithm for the Design of Tuned Mass Dampers under Earthquake Excitations. *Appl. Sci.* **2021**, *11*, 6172. https://doi.org/10.3390/app11136172

Academic Editor: Eui-Nam Huh

Received: 28 May 2021 Accepted: 28 June 2021 Published: 2 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

have been proposed thereafter, considering the harmonic loads, and random stationary and nonstationary white noise processes for single degree of freedom systems [6–9].

The first full-scale implementations of TMDs were aimed to control dynamic displacements caused by wind-induced vibrations. In that sense, TMDs have been deployed in several structures around the world, including the Centre Point Tower in Sydney, the Citicorp Center in New York, the Chiba Port Tower in Japan, the CN Tower in Canada, and, more recently, the Taipei 101 tower in Taiwan [10,11]. Several variations based on the conventional TMD formulation have been studied and applied thereafter, including pendulum-tuned mass dampers (PTMD) [12–14], tuned liquid column dampers (TLCD) [15–17], bidirectional TMDs [18] and hybrid, semi-active, and active TMDs [19–21]. However, some early investigations concluded that TMDs were not effective in reducing the response of buildings subjected to seismic excitation [22,23], which was due to two principal reasons. First, the limitations associated with the amount of mass added to the structural system, and, more importantly, that conventional TMDs can suppress effectively a single vibration mode. Instead, wide-band multi-modal damping can be achieved by complex systems, such as nonlinear TMDs [24,25], multiple positioned TMDs [26], or the amplification of the inertial mass in conventional TMDs using inerter devices (TMDIs) [27–29].

In spite of this, since the early 2000s, some researchers have proved the effectiveness of numerical iterative methods in the tuning of linear mass dampers for seismic applications. Classic techniques have been used for that purpose, such as minimax optimization [30], response surface methodology [31], and nonlinear gradient-based optimization [32]. Moreover, conventional population-based metaheuristics and evolutionary algorithms have been used to improve the computational efficiency in the tuning procedure; among them, particle swarm optimization [33,34]. harmony search [35–40], ant colony [41,42], flower pollination algorithm [43–45], bat algorithm [46], cuckoo search [47], genetic algorithms [48–50], and, more recently, chaotic optimization [51]. Furthermore, various of these investigations have included the effects of dynamic soil–structure interaction (SSI) [52–55]. Even though the tuning of linear TMDs on single degree of freedom systems is a very well addressed problem in the literature, most of these works were based on the assumption that closedform expressions, like those presented in [6–9], are not valid for multi-degree of freedom systems subjected to actual seismic loads. On the contrary, the usage of metaheuristic techniques or evolutionary algorithms allows the best-fit design variables for linear TMDs to be determined, using actual accelerograms as input excitations.

In that sense, this study introduces a novel methodology for the optimal design of passive TMDs located at upper levels of high-rise buildings to improve the seismic safety of structures based on the Whale Optimization Algorithm (WOA) [56]. The algorithm is modified to reduce the maximum horizontal peak displacement of the structure, and the root mean square (RMS) response of displacements as well. Moreover, four additional objective functions, derived from multiple weighted linear combinations of these parameters, are also studied in order to obtain the most efficient TMD design configuration. The results from the WOA optimization are compared with the differential evolution method (DEM) [57], whose effectiveness has been demonstrated extensively for TMDs and TMDIs seismic applications in previous works [58,59], and an exhaustive search (ES) process with precision to two decimal positions. Then, the proposed methodology is applied to a 32- story case-study derived from an actual building structure, and various accelerograms of recorder earthquakes are considered to assess its seismic performance in the linear-elastic range. The results of the study show a large enhancement in the dynamic response of the building controlled by the WOA designed TMDs, achieving reductions in the maximum floor displacements of up to 43%. Finally, an attempt to establish a single set of design parameters for TMDs based on the methodology proposed by Fallah and Zamiri [60] is implemented to verify the robustness of the optimization.

#### **2. Equations of Motion of** *n***-DOF Building Structures Controlled via TMDs for Seismic Applications**

The equations of motion that govern the behavior of the *n*-DOF structural system equipped with a TMD located at the top floor of the building are deduced from the scheme depicted in Figure 1. The displacement at each story-level is represented by *xi*, while the stiffness, damping, and mass coefficients of each story-level of the building are defined as *ki*, *ci*, and *mi*, respectively. Besides, the TMD action is introduced as an extra DOF represented by *xd*, while the stiffness, damping, and mass parameters are identified as *kd*, *cd*, and *md*, respectively.

Subsequently, the structural system controlled via TMD can be modeled as a set of *n* + 1 degrees of freedom, whose equation of motion can be expressed as:

$$\mathbf{M}\ddot{\mathbf{x}}(t) + \mathbf{C}\dot{\mathbf{x}}(t) + \mathbf{K}\mathbf{x}(t) = -\mathbf{M}1\ddot{\mathbf{x}}\_{\mathcal{K}}(t) \tag{1}$$

where, **x**(*t*), **ẋ**(*t*), and **ẍ**(*t*) are the *n* + 1-order vectors of displacement, velocity, and acceleration of the system; **1** is a unit vector, and *<sup>ẍ</sup><sup>g</sup>g*(*t*) denotes the ground acceleration over time. On the other hand, the matrices **M**, **C,** and **K** represent the mass, damping, and stiffness matrices of the system, which can be expressed as:

$$\mathbf{M} = \begin{bmatrix} m\_1 & 0 & 0 & 0 & \cdots & \cdots & 0 \\ 0 & m\_2 & 0 & 0 & \cdots & \cdots & 0 \\ 0 & 0 & m\_3 & 0 & \cdots & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \cdots & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & m\_{n-1} & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & 0 & m\_n & 0 \\ 0 & 0 & \cdots & \cdots & 0 & 0 & m\_d \end{bmatrix} \tag{2}$$


It should be noted that matrices (2) to (4) are idealized as tridiagonal stiffness matrices typical of shear frame structures. Nevertheless, the main purpose of these mathematical expressions is to illustrate the inclusion of the TMD effect on a linear structural system with *n* horizontal degrees of freedom, and, therefore, it can be extended to other structural models such as 2-dimensional frames derived from actual building numerical models with a static condensation applied to all other vertical and rotational degrees of freedom. Now, the equation of motion of the system is modified to a space–state representation to determine the dynamic response of the linear system. Hence, the state vector **z**(*t*) is introduced as: 

$$\mathbf{z}(t) = \left\{ \begin{array}{c} \mathbf{x}(t) \\ \dot{\mathbf{x}}(t) \end{array} \right\} \tag{5}$$

The above-mentioned vector involves the displacements and velocities of the controlled system. By deploying Equation (1), the state–space representation of the structural system controlled by TMD can be written as:

$$\left\{ \begin{array}{c} \dot{\mathbf{x}}(t) \\ \ddot{\mathbf{x}}(t) \end{array} \right\} = \left[ \begin{array}{c} 0 \\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{C} \end{array} \right] \left[ \begin{array}{c} \mathbf{x}(t) \\ \dot{\mathbf{x}}(t) \end{array} \right] + \left[ \begin{array}{c} 0 \\ \mathbf{M}^{-1} \end{array} \right] \cdot -\mathbf{M}1\ddot{\mathbf{x}}\_{\mathcal{S}}(t) \tag{6}$$

where **0** and **I** denote a zero matrix and an identity matrix, respectively. Subsequently, Equation (6) is rewritten as:

$$
\dot{\mathbf{z}}(t) = \mathbf{A}z(t) + \mathbf{B} \cdot - \mathbf{M}1\ddot{\mathbf{x}}\_{\mathcal{K}}(t) \tag{7}
$$

where **A** represents the transition state matrix, and **B** the location adjustment matrix of the external excitation in the structural system:

$$\mathbf{A} = \begin{bmatrix} 0 & \mathbf{I} \\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{C} \end{bmatrix} \tag{8}$$

$$\mathbf{B} = \begin{bmatrix} 0 \\ \mathbf{M}^{-1} \end{bmatrix} \tag{9}$$

Taking inspiration from the tuning process of TMDs on single degree of freedom systems, and to simplify the optimization process, the *kd* and *cd* parameters can be rewritten as:

$$k\_d = \omega\_s^2 f^2 m\_d \tag{10}$$

$$c\_d = 2\zeta\_d f \omega\_{\mathfrak{s}} m\_d \tag{11}$$

with *ω*s as the natural frequency of the structural system and *f* and *ζd* as the dimensionless frequency and damping ratios to be optimized. The reference value *ω*s will be assumed in this investigation as the corresponding frequency to the first vibration mode of the structure. Nevertheless, this is just an assumption to reduce significantly the search domain of the variables *kd* and *cd* since the optimization problem is focused on the tuning problem for TMDs on multi-degree-of-freedom systems.

#### **3. Whale Optimization Algorithm**

The Whale Optimization Algorithm (WOA) is a swarm-based metaheuristic recently proposed by Mirjalili and Lewis in [56], who developed the algorithm by taking inspiration from the bubble-net hunting strategy typically employed by humpback whales. In the last few years, other researchers have applied multiple modifications to the algorithm in order to solve various engineering optimization problems. For instance, Kaveh and Ghazaan [61] studied the sizing optimization of skeletal structures, Chen et al. [62] solved bar truss design and I-beam design problems, and Azizi et al. [63] optimized a fuzzy controller for seismically excited nonlinear structures.

WOA strategy first defines a population of whales and evaluates a set of random solutions according to both the position of each whale and a certain objective function. By comparing the primary random solutions, WOA selects an initial best search agent. In each iteration, whales update their position with respect to either a randomly chosen search agen<sup>t</sup> or the best solution obtained so far. Then, the new positions are calculated according to Equations (12) and (13):

$$\Delta = |\Gamma \cdot \mathbf{X}^\*(w) - \mathbf{X}(w)|\tag{12}$$

$$\mathbf{X}(w+1) = \mathbf{X}^\*(w) - \mathbf{Y} \cdot \boldsymbol{\Delta} \tag{13}$$

where *w* is the current whale generation, and **X**\* is the updated position vector of the best solution. In addition, the coefficient vectors Ψ and Γ are calculated as follows:

$$\mathbf{Y} = 2a \cdot \mathbf{r} - a \tag{14}$$

$$
\Gamma = 2\mathbf{r} \tag{15}
$$

with **r** as a random value in the domain (0, 1). Besides, *a* denotes a sequentially decreasing number from 2 to 0 as whales generation *w* increases, to provide exploration and exploitation.

According to Mirjalili and Lewis [56], a random search agen<sup>t</sup> is chosen when | Ψ| ≥ 1, while the best solution is selected when | Ψ| < 1 for updating the position of the search agents. The bubble-net feeding behavior of humpback whales serves as an inspiration to update the position of the search agents. It is assumed that there is the same probability of selecting between either a shrinking encircling mechanism or a spiral model, to improve the position of the search agents during the optimization procedure and ge<sup>t</sup> closer to the optimal solution. This concept is mathematically formulated as:

$$\mathsf{X}(w+1) = \begin{cases} \mathsf{X}^\*(w) - \mathsf{Y} \cdot \Delta & \text{if } p < 0.5\\ \Delta^\prime \mathsf{e}^{bl} \cdot \cos(2\pi l) + \mathsf{X}^\*(w) & \text{if } p \ge 0.5 \end{cases} \tag{16}$$

where *b* is a constant that defines the shape of the logarithmic spiral; *l* is a random number in the space ( −1, 1); *p* is a random probability in (0, 1). In the adaptation of the algorithm, the constant *b* was set to 1 following the recommendations provided in [61–63]. On the other hand, the *l* and *p* parameters are chosen randomly in their respective domain for every generation *w*, in order to avoid convergence at local minimums and to improve the efficiency of the algorithm. Finally, Δ defines the distance of the *i*th whale to the best solution obtained so far, which could be calculated as:

$$
\Delta' = |\mathbf{X}^\*(w) - \mathbf{X}(w)|\tag{17}
$$

Six objective functions are proposed in the optimization process to diminish the dynamic response of the structure when it is subjected to seismic excitations. These functions are associated with the reduction of the maximum horizontal peak of displacement of the structure, and the reduction of the root mean square (RMS) values of displacements. Certainly, there are other critical parameters to be analyzed in a robust structural design, such as inter-story drifts and peak floor accelerations, which are closely related to structural damage. Nevertheless, the criteria for selecting the horizontal peak displacements and RMS values of displacements as critical parameters in the tuning procedure responds to two principal reasons. Firstly, the reduction of absolute displacements in the objective functions contributes directly to diminishing the inter-story drift values and, more importantly, the main objective of this research is to propose a novel tuning methodology based on WOA, as well as proving its efficiency over other optimization techniques [30–55]. Moreover, the usage of multiple functions aims to determine which of these objective functions is more efficient, based on the evaluation and comparison on the decrease in the response of both previously mentioned parameters.

In that sense, the objective functions OA1 and OA2 correspond to the reduction of the maximum horizontal peak displacement in the entire structure, and the decrease in the RMS values of displacements, respectively:

$$\text{OA1} = \min(\max(|\mathbf{x}\_n|)) \quad \text{for } \mathbf{x}\_n = [\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n] \tag{18}$$

$$\text{OA2} = \min(\max(RMS(\mathbf{x}\_n))) \quad \text{for } \mathbf{x}\_n = [\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n] \tag{19}$$

Regarding the remaining four objective functions, these are defined as a variable weighted linear combination of the relationships between the above-mentioned parameters in Equations (18) and (19). These functions are defined as J1, J2, J3, and J4, and are described through expression (20):

$$\mathbf{J}\_{i} = \min\left( (i \ast 0.2) \frac{\max(|\mathbf{x}\_{n}|)}{\max\left(|\mathbf{x}\_{n}^{\*}|\right)} + (1 - (i \ast 0.2)) \frac{\max(\text{RMS}(\mathbf{x}\_{n}))}{\max\left(\text{RMS}(\mathbf{x}\_{n}^{\*})\right)} \right) \text{ for } i = 1:4 \text{ and } \mathbf{x}\_{n} = [\mathbf{x}\_{1}, \ \mathbf{x}\_{2}, \dots, \mathbf{x}\_{n}] \tag{20}$$

where x<sup>∗</sup> n. denotes the uncontrolled response of the structure at the *n* degree of freedom of the structure, introduced to normalize the values of J*i* function. Different from OA1 and OA2 functions, in which absolute values can be used to analyze the numerical results computed from the optimization process, the function J*i* combines horizontal peak displacements and RMS values of displacements which are commonly different in magnitude. Hence, both values must be normalized to be properly weighted and combined thereafter.

The TMD tuning procedure in multi-degree of freedom systems is a two-dimensional optimization problem focused on finding the optimal *kd* and *cd* values for a fixed mass ratio (*μ*). As mentioned in Section 2, the optimization process is simplified by the inclusion of Equations (10) and (11). Hence, WOA is programed to find the optimal *f* and *ζd* values in the domain:

$$0.50 \le f \le 2.00\tag{21}$$

$$0 \le \mathcal{Z}\_d \le 0.50 \tag{22}$$

The bounds in the search domain were defined according to practical recommendations that were taken from [30–55]. The flowchart in Figure 2 summarizes the WOA procedure.

**Figure 2.** Description of the WOA procedure.
