A Novel Hybrid Harmony Search Approach for the Analysis of Plane Stress Systems via Total Potential Optimization

Abstract

By finding the minimum total potential energy of a structural system with a defined degree of freedoms assigned as design variables, it is possible to find the equilibrium condition of the deformed system. This method, called total potential optimization using metaheuristic algorithms (TPO/MA), has been verified on truss and truss-like structures, such as cable systems and tensegric structures. Harmony Search (HS) algorithm methods perfectly found the analysis results of the previous structure types. In this study, TPO/MA is presented for analysis of plates for plane stress members to solve general types of problems. Due to the complex nature of the system, a novel hybrid Harmony Search (HHS) approach was proposed. HHS is the hybridization of local search phases of HS and the global search phase of the Flower Pollination Algorithm (FPA). The results found via HHS were verified with the finite element method (FEM). When compared with classical HS, HHS provides smaller total potential energy values, and needs less iterations than other new generation metaheuristic algorithms.

1. Introduction

In the analysis of structural members, the general trend for linear systems is to find deflections after forming a matrix equation involving the stiffness matrix, loads, and deflections. By using very elaborate matrix methods, linear equations thus formed are solved with great success. This is not the case for nonlinear systems, where the use of additional iterative methods becomes a must [1,2,3,4] since the stiffness matrix halts to be a constant one. Effectively, in those cases, the stiffness function is in turn a function of displacements. Detailed information about the optimization and design of various structure types can be found in Nazir et al. [5].

In recent years, total potential optimization using metaheuristic algorithms (TPO/MA) has been proposed as an iterative analysis method [6]. As is already known, all mechanical systems are in stable equilibrium in the case of minimum total potential energy. Thus, the minimization of the total potential energy of a system can be performed by considering the deflections as design variables. TPO/MA have been used successfully for the analysis of several structural systems. The first application area of TPO/MA is the truss structures. Different metaheuristic algorithms, such as local search [6], harmony search (HS) [7], and particle swarm optimization (PSO) [8], are employed with TPO/MA to analyze truss structures.

In TPO/MA, since the displacements are the basic unknowns of the problem, the method becomes effective in finding second-order effects directly without the requirement of an additional process. Additionally, the finding of non-unique solutions [9] and material non-linearity [10] can be easily tackled.

Different to truss structures, other truss-like systems, such as tensegric structures [11] and cable systems [12], are also solved via TPO/MA. In order to enhance the capacity of the method in solving large scale systems, hybrid metaheuristic algorithms combining the flower pollination algorithm (FPA), the teaching learning-based optimization (TLBO), and the Jaya algorithm (JA) are also proposed [13].

In this study, an innovation has been performed for enlarging the application area of the method TPO/MA to plates and plate-like structures. For this purpose, a simple approach has been adopted, which consists of dividing the system into triangles, of which the nodal displacements are the unknowns of the problem. Then, the strain energy of each triangle is expressed in terms of these variables. The summation of all these energies yields the total strain energy of the system. Adding to this value algebraically, the work done by the external forces again expressed in terms of the nodal deflections yield the total potential energy of the whole system considered. The rest of the technique is the determination of nodal displacements that will minimize the total potential energy through a metaheuristic optimization algorithm. Only plane stress problems are considered in this first trial.

According to previous studies, HS is a suitable and effective algorithm for TPO/MA. For the numerical examples of the study, HS is effective to find close values to the finite element method (FEM) results and other new generated algorithms. Additionally, HS is also a well-proved algorithm for structural optimization problems [14,15]. Several recent applications in the last decade are mentioned as follows.

Kaveh and Talatahari developed a hybrid algorithm that is based on PSO with passive congregation, ant colony optimization, and HS to solve the optimum design of plane and space truss structures [16]. An adaptive HS algorithm was proposed by Hasançebi et al., and the proposed method for the optimum design of large-scale steel frame structures has the ability to automatically adjust the algorithm parameters [17]. Martini proposed HS for structural frameworks for size, shape and topology optimization in order to conduct a multimodal optimization proposing a range of good alternative solutions instead of a single optimum one [18]. Degertekin proposed two different modifications of HS for the optimization of truss structures [19]. The sequential HS algorithm was proposed to find the shape optimization of structural systems, including frequency constraints [20]. For the control of the structures, HS-based methodologies are also perfectly applied for the placement of dampers [21], tuned mass dampers [22,23], and base isolators [24]. In structural health monitoring, Jin et al. proposed and improved HS for optimal sensor placement to crane structures [25]. For the optimization of double- and triple-layer grid structures, a modified gravitational HS algorithm was used with the efficient member grouping strategy [26]. Elaki et al. proposed a method combining HS and Artificial Neural Networks for the aim of assessing the reliability of structures with the implicit limit state function [27]. By using the three-parameter Vlaslov Elastic Foundation model, soil–structure interaction was also considered in the optimum design of braced steel frame structures using HS and TLBO [28]. De Almedia investigated the optimum design of laminated composite structures via HS [29]. The modification of HS for the structure of the algorithm is presented in a review study conducted by Zhang and Geem [30].

The main criticized factor of TPO/MA is the use of several iterations, while FEM can find the results directly for linear systems. Although FEM also needs the iteration methods for non-linear systems, the number of iterations is the most important doubt for TPO/MA. The duration of analysis is generally short for TPO/MA, but the area of using algorithms with a fast convergence ability is still a needed and important area. Due to that, hybrid algorithms are needed to develop for sensitive and fast analysis. A novel hybrid approach was proposed to find comparative results that can be found with less iterations than other algorithms. New generation metaheuristic algorithms are easy to apply to engineering problems directly and the number of parameters is generally low, and these parameters are not critically effective on solutions. Due to that, Harmony Search (HS) and the Flower Pollination Algorithm (FPA), with a Lévy distribution different from a classical random generation between 0 and 1, are hybridized instead of using classical algorithms, such as Genetic Algorithms and PSO. The hybrid Harmony Search (HHS) approach combined with FPA is an effective tool for this goal according to numerical experiments proposed in the study.

2. The Total Potential Energy of Plane Stress Members

As explained in the above paragraphs, the plates considered are divided into triangular elements such that the state variables are the displacements u_i and v_i at the nodes in x and y directions, respectively. Considering a linear variation in the displacements over the triangle, the displacement fields, such as u(x,y) and v(x,y), can be written as in Equations (1) and (2), respectively.

$$\mathrm{u}(\mathrm{x},\mathrm{y})={\mathrm{u}}_{\mathrm{i}}+{\mathrm{C}}_1\mathrm{x}+{\mathrm{C}}_2\mathrm{y},$$

(1)

$$\mathrm{v}(\mathrm{x},\mathrm{y})={\mathrm{v}}_{\mathrm{i}}+{\mathrm{C}}_3\mathrm{x}+{\mathrm{C}}_4\mathrm{y},$$

(2)

When the particular derivatives of displacement fields are taken with respect to x and y, the strains are expressed in terms of the constants C₁, C₂, C₃, and C₄.

$${\epsilon }_{\mathrm{x}}=\frac{\partial \mathrm{u}}{\partial \mathrm{x}}={\mathrm{C}}_1,$$

(3)

$${\epsilon }_{\mathrm{y}}=\frac{\partial \mathrm{v}}{\partial \mathrm{y}}={\mathrm{C}}_4,$$

(4)

$$\gamma =\frac{\partial \mathrm{u}}{\partial \mathrm{y}}+\frac{\partial \mathrm{v}}{\partial \mathrm{x}}={\mathrm{C}}_2+{\mathrm{C}}_3,$$

(5)

The strains ε_x, ε_y, and γ define the normal strain in x and y directions and the shear strain, respectively. If the coordinates of node i are the origin point, the nodal displacements of nodes i, j, and k can be expressed as follows, according to Figure 1:

$$\mathrm{u}(0,0)={\mathrm{u}}_{\mathrm{i}},\hspace{0.17em}\hspace{0.17em}\mathrm{v}(0,0)={\mathrm{v}}_{\mathrm{i}},$$

(6)

$$\mathrm{u}({\mathrm{a}}_{\mathrm{j}},{\mathrm{b}}_{\mathrm{j}})={\mathrm{u}}_{\mathrm{j}},\hspace{0.17em}\mathrm{v}({\mathrm{a}}_{\mathrm{j}},{\mathrm{b}}_{\mathrm{j}})={\mathrm{v}}_{\mathrm{j}},$$

(7)

$$\mathrm{u}({\mathrm{a}}_{\mathrm{k}},{\mathrm{b}}_{\mathrm{k}})={\mathrm{u}}_{\mathrm{k}},\hspace{0.17em}\mathrm{v}({\mathrm{a}}_{\mathrm{k}},{\mathrm{b}}_{\mathrm{k}})={\mathrm{v}}_{\mathrm{k}},$$

(8)

By using Equations (1) and (2), the relations between the nodal displacements is written as Equation (9):

$$\left[\begin{array}{c}{\mathrm{u}}_{\mathrm{j}}\\ {\mathrm{v}}_{\mathrm{j}}\\ {\mathrm{u}}_{\mathrm{k}}\\ {\mathrm{v}}_{\mathrm{k}}\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{a}}_{\mathrm{j}}& {\mathrm{b}}_{\mathrm{j}}& 0& 0\\ 0& 0& {\mathrm{a}}_{\mathrm{j}}& {\mathrm{b}}_{\mathrm{j}}\\ {\mathrm{a}}_{\mathrm{k}}& {\mathrm{b}}_{\mathrm{k}}& 0& 0\\ 0& 0& {\mathrm{a}}_{\mathrm{k}}& {\mathrm{b}}_{\mathrm{k}}\end{array}\right]\left[\begin{array}{c}{\mathrm{C}}_1\\ {\mathrm{C}}_2\\ {\mathrm{C}}_3\\ {\mathrm{C}}_4\end{array}\right]+\left[\begin{array}{c}{\mathrm{u}}_{\mathrm{i}}\\ {\mathrm{v}}_{\mathrm{i}}\\ {\mathrm{u}}_{\mathrm{i}}\\ {\mathrm{v}}_{\mathrm{i}}\end{array}\right],$$

(9)

The solution of Equation (9) yields the constants C₁–C₄ as in Equation (10)–(13):

$${\mathrm{C}}_1=\frac{{\mathrm{b}}_{\mathrm{k}}({\mathrm{u}}_{\mathrm{j}}-{\mathrm{u}}_{\mathrm{i}})}{{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}-{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}}+\frac{{\mathrm{b}}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{k}}-{\mathrm{u}}_{\mathrm{i}})}{{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}-{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}},$$

(10)

$${\mathrm{C}}_2=\frac{{\mathrm{a}}_{\mathrm{k}}({\mathrm{u}}_{\mathrm{j}}-{\mathrm{u}}_{\mathrm{i}})}{{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}-{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}}+\frac{{\mathrm{a}}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{k}}-{\mathrm{u}}_{\mathrm{i}})}{{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}-{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}},$$

(11)

$${\mathrm{C}}_3=\frac{{\mathrm{b}}_{\mathrm{k}}({\mathrm{v}}_{\mathrm{j}}-{\mathrm{v}}_{\mathrm{i}})}{{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}-{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}}+\frac{{\mathrm{b}}_{\mathrm{j}}({\mathrm{v}}_{\mathrm{k}}-{\mathrm{v}}_{\mathrm{i}})}{{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}-{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}},$$

(12)

$${\mathrm{C}}_4=\frac{{\mathrm{a}}_{\mathrm{k}}({\mathrm{v}}_{\mathrm{j}}-{\mathrm{v}}_{\mathrm{i}})}{{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}-{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}}+\frac{{\mathrm{a}}_{\mathrm{j}}({\mathrm{v}}_{\mathrm{k}}-{\mathrm{v}}_{\mathrm{i}})}{{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}-{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}},$$

(13)

The strain energy density that resulted due to deformation is written as Equation (14) for a two-dimensional elastic body:

$$\mathrm{e}={\displaystyle \underset{\epsilon =0}{\overset{\epsilon }{\mathrm{\int }}}\mathsf{\sigma }\mathrm{d}\epsilon =\frac{1}{2}({\mathsf{\sigma }}_{\mathrm{x}}{\epsilon }_{\mathrm{x}}}+{\mathsf{\sigma }}_{\mathrm{y}}{\epsilon }_{\mathrm{y}}+\mathsf{\tau }\mathsf{\gamma }),$$

(14)

When the solution of constants given as Equations (10)–(13) are used in Equations (3)–(5), the stains are written as follows:

$${\epsilon }_{\mathrm{x}}=\frac{{\mathrm{b}}_{\mathrm{k}}({\mathrm{u}}_{\mathrm{j}}-{\mathrm{u}}_{\mathrm{i}})}{{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}-{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}}+\frac{{\mathrm{b}}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{k}}-{\mathrm{u}}_{\mathrm{i}})}{{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}-{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}},$$

(15)

$${\epsilon }_{\mathrm{y}}=\frac{{\mathrm{a}}_{\mathrm{k}}({\mathrm{v}}_{\mathrm{j}}-{\mathrm{v}}_{\mathrm{i}})}{{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}-{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}}+\frac{{\mathrm{a}}_{\mathrm{j}}({\mathrm{v}}_{\mathrm{k}}-{\mathrm{v}}_{\mathrm{i}})}{{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}-{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}},$$

(16)

$$\mathsf{\gamma }=\frac{1}{{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}-{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}}+\left({\mathrm{a}}_{\mathrm{k}}({\mathrm{u}}_{\mathrm{j}}-{\mathrm{u}}_{\mathrm{i}})+a_{\mathrm{j}}({\mathrm{u}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{k}})+{\mathrm{b}}_{\mathrm{k}}({\mathrm{v}}_{\mathrm{i}}-{\mathrm{v}}_{\mathrm{j}})+{\mathrm{b}}_{\mathrm{j}}({\mathrm{v}}_{\mathrm{k}}-{\mathrm{v}}_{\mathrm{i}})\right),$$

(17)

According to linear stress–strain relations for two-dimensional plane stress problems, stresses can be found from Equations (18)–(20), where E is the elasticity modulus and ν is the Poisson’s ratio:

$${\mathsf{\sigma }}_{\mathrm{x}}=\frac{\mathrm{E}}{1-{\nu }^2}({\epsilon }_{\mathrm{x}}+\nu {\epsilon }_{\mathrm{y}}),$$

(18)

$${\mathsf{\sigma }}_{\mathrm{y}}=\frac{\mathrm{E}}{1-{\nu }^2}({\epsilon }_{\mathrm{y}}+\nu {\epsilon }_{\mathrm{x}}),$$

(19)

$$\mathsf{\tau }=\frac{\mathrm{E}}{1-{\nu }^2}(\frac{1-\nu }{2}\mathsf{\gamma }),$$

(20)

The substitution of Equations (15)–(17) to (18)–(20) yields Equations (21)–(23) for stresses for linear problems:

$${\mathsf{\sigma }}_{\mathrm{x}}=\frac{\mathrm{E}}{({\nu }^2-1)({\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}-{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}})}\nu \left({\mathrm{a}}_{\mathrm{k}}({\mathrm{v}}_{\mathrm{i}}-{\mathrm{v}}_{\mathrm{j}})+{\mathrm{a}}_{\mathrm{j}}({\mathrm{v}}_{\mathrm{k}}-{\mathrm{v}}_{\mathrm{i}})\right)+(\nu -1)\left({\mathrm{b}}_{\mathrm{k}}({\mathrm{u}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{j}})-{\mathrm{b}}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{k}})\right),$$

(21)

$${\mathsf{\sigma }}_{\mathrm{y}}=\frac{\mathrm{E}}{({\nu }^2-1)({\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}-{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}})}{\mathrm{a}}_{\mathrm{k}}({\mathrm{v}}_{\mathrm{i}}-{\mathrm{v}}_{\mathrm{j}})+{\mathrm{a}}_{\mathrm{j}}({\mathrm{v}}_{\mathrm{k}}-{\mathrm{v}}_{\mathrm{i}})+\nu {\mathrm{b}}_{\mathrm{k}}({\mathrm{u}}_{\mathrm{j}}-{\mathrm{u}}_{\mathrm{i}})+\nu {\mathrm{b}}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{k}}),$$

(22)

$$\mathsf{\tau }=\frac{\mathrm{E}}{({\nu }^2-1)({\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}-{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}})}(\frac{1-\nu }{2}{\mathrm{a}}_{\mathrm{k}}({\mathrm{u}}_{\mathrm{j}}-{\mathrm{u}}_{\mathrm{i}})+{\mathrm{a}}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{k}})+{\mathrm{b}}_{\mathrm{k}}({\mathrm{v}}_{\mathrm{i}}-{\mathrm{v}}_{\mathrm{j}})+{\mathrm{b}}_{\mathrm{j}}({\mathrm{v}}_{\mathrm{k}}-{\mathrm{v}}_{\mathrm{i}}),$$

(23)

Then, the strain energy for linear problems can be found from Equation (14), as can be seen in Equation (24), for one element:

$$\begin{array}{l}\mathrm{e}=\frac{\mathrm{E}}{4(1-{\nu }^2){({\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}}-{\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}})}^2}(-(\nu -1\left)\right(({\mathrm{a}}_{\mathrm{k}}({\mathrm{u}}_{\mathrm{j}}-{\mathrm{u}}_{\mathrm{i}})+{\mathrm{a}}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{k}})+{\mathrm{b}}_{\mathrm{k}}({\mathrm{v}}_{\mathrm{i}}-{\mathrm{v}}_{\mathrm{j}})+{\mathrm{b}}_{\mathrm{j}}({\mathrm{v}}_{\mathrm{k}}-{\mathrm{v}}_{\mathrm{i}}){)}^2\\ +2({\mathrm{a}}_{\mathrm{k}}({\mathrm{v}}_{\mathrm{j}}-{\mathrm{v}}_{\mathrm{i}})+{\mathrm{a}}_{\mathrm{j}}({\mathrm{v}}_{\mathrm{i}}-{\mathrm{v}}_{\mathrm{k}}))(\nu ({\mathrm{a}}_{\mathrm{k}}({\mathrm{v}}_{\mathrm{j}}-{\mathrm{v}}_{\mathrm{i}})+{\mathrm{a}}_{\mathrm{j}}({\mathrm{v}}_{\mathrm{i}}-{\mathrm{v}}_{\mathrm{k}}))+{\mathrm{b}}_{\mathrm{k}}({\mathrm{u}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{j}})+{\mathrm{b}}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{k}}-{\mathrm{u}}_{\mathrm{i}}))\\ +2({\mathrm{a}}_{\mathrm{k}}({\mathrm{v}}_{\mathrm{j}}-{\mathrm{v}}_{\mathrm{i}})+{\mathrm{a}}_{\mathrm{j}}({\mathrm{v}}_{\mathrm{i}}-{\mathrm{v}}_{\mathrm{k}}))({\mathrm{a}}_{\mathrm{k}}({\mathrm{v}}_{\mathrm{j}}-{\mathrm{v}}_{\mathrm{i}})+{\mathrm{a}}_{\mathrm{j}}({\mathrm{v}}_{\mathrm{i}}-{\mathrm{v}}_{\mathrm{k}}))+\nu {\mathrm{b}}_{\mathrm{k}}({\mathrm{u}}_{\mathrm{i}}-{\mathrm{u}}_{\mathrm{j}})+\nu {\mathrm{b}}_{\mathrm{j}}({\mathrm{u}}_{\mathrm{k}}-{\mathrm{u}}_{\mathrm{i}})))\end{array}$$

(24)

The strain energy density of the m^th element (e_m) is multiplied with the volume (V_m) of the element to yield the strain energy of the m^th element, as in Equation (25):

$${\mathrm{U}}_{\mathrm{i}}={\mathrm{e}}_{\mathrm{m}}{\mathrm{V}}_{\mathrm{m}},$$

(25)

$${\mathrm{V}}_{\mathrm{m}}=\frac{({\mathrm{a}}_{\mathrm{j}}{\mathrm{b}}_{\mathrm{k}}-{\mathrm{a}}_{\mathrm{k}}{\mathrm{b}}_{\mathrm{j}})\mathrm{t}}{2},$$

(26)

The volume is the multiplication of thickness (t) and area of the triangle element, as shown in Equation (26). The total strain energy of the system is calculated by adding the strain energy of all the elements. For the system consisting of n elements, the strain energy is as follows:

$$\mathrm{U}={\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{n}}{\mathrm{U}}_{\mathrm{m}}},$$

(27)

In order to find the total potential energy (Π_p), the work done by the external forces must be subtracted from the total strain energy. Equation (28) is given for a system with p number of nodes and point loads in x and y directions (P_xn and P_yn):

$${\mathsf{\Pi }}_{\mathrm{p}}=\mathrm{U}-{\displaystyle \sum _{\mathrm{n}=1}^{\mathrm{p}}({\mathrm{P}}_{\mathrm{x}\mathrm{n}}{\mathrm{u}}_{\mathrm{n}}+{\mathrm{P}}_{\mathrm{y}\mathrm{n}}{\mathrm{v}}_{\mathrm{n}})},$$

(28)

3. The Analysis Methodology

In TPO/MA, the nodal displacements (u_n and v_n for m=1 to p) are assigned with candidate values according to the rules of employed metaheuristic algorithms. Each node has two nodal displacements and the number of design variables (s) is equal to two times the number of nodes (p) minus the number of fixed displacements (p_f):

$$\mathrm{s}=2\mathrm{p}-{\mathrm{p}}_{\mathrm{f}},$$

(29)

The flowchart of the TPO/MA method for plane stress members is shown as Figure 2. The optimization process starts with the definition of the problem. Then, the nodal displacements must be assumed to find the total potential energy of the system. The initial solution matrix is generated via the candidate solutions of nodal displacements, and the assumed values are randomly chosen from the user-defined range of nodal displacements. Then, the iterative and essential optimization process starts, and special algorithm rules are taken into consideration for the maximum number of iterations.

In the methodology, there are two phases. Generally, one of these phases is global optimization, while the other one is a local optimization. These two phases are chosen via a parameter, and it is compared with a randomly generated number between 0 and 1 (R). The two phases of HS and HHS are chosen according to a parameter called the harmony memory considering rate (HMCR).

In this paper, the music inspired Harmony Search (HS) algorithm developed by Geem et al. [31] was employed and it was hybridized with the Flower Pollination Algorithm (FPA) developed by Yang [32]. The results of numerical examples were also compared with single forms of HS and the FPA, the Jaya Algorithm (JA) [33], and teaching–learning-based optimization (TLBO) [34].

After each iteration, the total potential energy values are calculated, and these values are compared with the existing values for updating.

3.1. Classical Harmony Search

In the classical form of HS, which has been applied to various structural optimization problems [35], the equation of the first phase is taken as Equation (30). This phase is the global search phase, which is also used for the generation of initial values.

$${\mathrm{X}}_{\mathrm{n}\mathrm{e}\mathrm{w}}={\mathrm{X}}_{\min }+\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}(0,1)({\mathrm{X}}_{\max }-{\mathrm{X}}_{\min }),$$

(30)

In Equation (30), the maximum and minimum ranges of design variables are shown as X_max and X_min, respectively. A random number between 0 and 1 (ran (0,1)) is generated to form a new set of design variables (X_new).

In local optimization phase of HS, a randomly selected existing solution (X_ran) is used, and the neighborhood range (X_ran ± fw) is scanned. The formulation of the local optimization of HS is shown as Equation (31):

$${\mathrm{X}}_{\mathrm{n}\mathrm{e}\mathrm{w}}={\mathrm{X}}_{\mathrm{r}\mathrm{a}\mathrm{n}}+\mathrm{f}\mathrm{w}(2\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}(0,1)-1)\hspace{0.17em}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\hspace{0.17em}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{a}\mathrm{b}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}\hspace{0.17em}\mathrm{o}\mathrm{f}\hspace{0.17em}\mathrm{P}\mathrm{A}\mathrm{R}$$

(31)

The parameter “pitch adjusting rate (PAR)” is a parameter controlling the acceptable range around the chosen solution.

3.2. Hybridization of Harmony Search

To increase the convergence ability of the classical HS approach, the best existing solution can be used in global optimization. Moreover, using a distribution method that is not linear distribution may be effective. For that reason, a Lévy distribution can be used instead of considering a random number between 0 and 1. The Flower Pollination Algorithm (FPA) developed by Yang uses Lévy distribution and the best solution in the global pollination phase, which is generated via the observations biotic and cross pollinations. Due to the behavior of pollinators, such as insects, birds and other animals, a random flight (Lévy flight) is considered. The formulation of the global optimization phase of FPA and Hybrid Harmony Search (HHS) is given as Equation (32). According to the formulation, the existing solution (X_old) is modified according to best exiting solution (g^*) and a Lévy distribution (L). Lévy distribution can be also assigned with random numbers covering a higher range than a linear distribution between 0 and 1. For that reason, it may affect the convergence ability negatively, but it will prevent to track to a local optimum result which is a close value to currently existing best solution that is not the global optimum.

$${\mathrm{X}}_{\mathrm{n}\mathrm{e}\mathrm{w}}={\mathrm{X}}_{\mathrm{o}\mathrm{l}\mathrm{d}}+\mathrm{L}({\mathrm{g}}^{*}-{\mathrm{X}}_{\mathrm{o}\mathrm{l}\mathrm{d}}),$$

(32)

4. Numerical Examples

As numerical examples, three problems were solved, and all solutions were compared with the finite element method (FEM) results. For all problems, the Harmony memory considering rate (HMCR) was taken as 0.5, while PAR is equal to 0.1. These parameters are the best cases for the examples in the study for sensitive final solutions and fast convergence. The HMCR value as 0.5 leads to usage of global and local optimization phases with equal probability, and the PAR value as 0.1 is a suitable value for convergence around existing solutions.

4.1. Structure 1: Cantilever Beam Under a Concentered Load

The figure of the first structure, which was previously considered by Cook [36,37], is given as Figure 3. The cantilever beam is loaded at the end (point C) in +y direction and the intensity of the load P is 20 N. The elasticity modulus and Poisson ratio of the material of the beam are 100 N/mm² and 0.3, respectively. The thickness of the members is 10 mm. The triangulation prepared for the analyses can be seen in the Figure 4.

The cantilever beam has 15 free nodes and 30 displacement values considered as design variables. As seen in the results presented in

Table 1. Solutions for the cantilever beam with 20 plane members

	TPOMA (FPA)		TPOMA (JA)		TPOMA (HS)		TPOMA (HHS)		TPOMA (TLBO)		FEM
Node	Δx (mm)	Δy (mm)	Δx (mm)	Δy (mm)	Δx (mm)	Δy (mm)	Δx (mm)	Δy (mm)	Δx (mm)	Δy (mm)	Δx (mm)	Δy (mm)
4	0.19369	0.25151	0.19369	0.25151	0.20715	0.27089	0.19259	0.24992	0.19369	0.25151	0.1937	0.2515
5	0.00059	0.24176	0.00059	0.24176	−0.00003	0.26078	0.00059	0.24022	0.00059	0.24176	0.0006	0.2418
6	−0.19912	0.25858	−0.19912	0.25858	−0.21326	0.27713	−0.19799	0.25693	−0.19912	0.25858	−0.1991	0.2586
7	0.34701	0.84814	0.34701	0.84814	0.36992	0.90709	0.34517	0.84317	0.34701	0.84814	0.347	0.8481
8	−0.00314	0.8425	−0.00314	0.8425	−0.0046	0.90067	−0.00312	0.83754	−0.00314	0.8425	−0.0031	0.8425
9	−0.35408	0.85627	−0.35408	0.85627	−0.37917	0.91413	−0.35220	0.85123	−0.35408	0.85627	−0.3541	0.8563
10	0.4547	1.70967	0.4547	1.70967	0.48581	1.83084	0.45247	1.70025	0.4547	1.70967	0.4547	1.7097
11	−0.00576	1.70704	−0.00576	1.70704	−0.00686	1.82991	−0.00571	1.69762	−0.00576	1.70704	−0.0058	1.707
12	−0.46619	1.71762	−0.46619	1.71762	−0.50155	1.8422	−0.46388	1.70816	−0.46619	1.71762	−0.4662	1.7176
13	0.51858	2.74696	0.51858	2.74696	0.55003	2.93925	0.51618	2.73258	0.51858	2.74696	0.5186	2.747
14	−0.0084	2.74754	−0.0084	2.74754	−0.01011	2.94063	−0.00834	2.73315	−0.0084	2.74754	−0.0084	2.7475
15	−0.53422	2.75474	−0.53422	2.75473	−0.57002	2.94848	−0.53172	2.74033	−0.53422	2.75473	−0.5342	2.7547
16	0.53861	3.87428	0.53861	3.87428	0.57192	4.13881	0.53619	3.85487	0.53861	3.87428	0.5386	3.8743
17	−0.01054	3.87495	−0.01054	3.87495	−0.01274	4.14002	−0.01048	3.85555	−0.01054	3.87495	−0.0105	3.8749
18	−0.55939	3.88219	−0.55939	3.88219	−0.59412	4.14664	−0.55686	3.86278	−0.55939	3.88219	−0.5594	3.8822
Пp (Nmm)	−38.8		−38.8		−38.6		−38.8		−38.8		−38.8
Analysis number	999103		672375		999343		531810		987341		−

, TPO/MA using HHS, FPA, JA, and TLBO has similar solutions to FEM. Application with HS resulted in a total potential energy value close to the others, but the nodal displacements have important differences, as compared to other algorithms. On the other hand, HS needs more iterations than other methods; this makes it computationally ineffective for the problem.

The single-phase JA has a better convergence rate as a classical form algorithm since it reached at minimum energy at 672375^th iteration, while the others were found nearly at the end of 10⁶ iterations. HHS has a significantly better convergence ability than JA by finding the best solution with 531810 iterations. By using HSS, all failures of HS are prevented.

4.2. Structure 2. A Plate with the Hole in the Middle (32 Members-25 Nodes)

The second structure is the plate with the hole in the middle [37,38], having a thickness of 2.5 mm (Figure 5). The material is characterized by the elasticity modulus 70000 N/mm² and Poisson ratio 0.3. The uniformly distributed load (q) with a density of 175 N/mm is applied to the plate in x direction. Because of symmetry, the quarter system shown in Figure 6 with 32 members and 25 nodes was solved.

The number of design variables is 40 because of the 10 fixed displacements of the 50 nodal displacements. The solutions of the design variables are shown in

Table 2. Solutions for the plate with the hole in the middle with 32 members.

	TPOMA (FPA)		TPOMA (JA)		TPOMA (HS)		TPOMA (HHS)		TPOMA (TLBO)		FEM
Node	Δx (mm)	Δy (mm)	Δx (mm)	Δy (mm)	Δx (mm)	Δy (mm)	Δx (mm)	Δy (mm)	Δx (mm)	Δy (mm)	Δx (mm)	Δy (mm)
1	0.000000	−0.006729	0.000000	−0.006729	0.000000	−0.006268	0.000000	−0.006729	0.000000	−0.006729	0.000000	−0.006729
2	0.000000	−0.010509	0.000000	−0.010509	0.000000	−0.010043	0.000000	−0.010509	0.000000	−0.010509	0.000000	−0.010509
3	0.000000	−0.014899	0.000000	−0.014899	0.000000	−0.014802	0.000000	−0.014899	0.000000	−0.014899	0.000000	−0.014899
4	0.000000	−0.021782	0.000000	−0.021782	0.000000	−0.021799	0.000000	−0.021782	0.000000	−0.021782	0.000000	−0.021782
5	0.000000	−0.033078	0.000000	−0.033078	0.000000	−0.032940	0.000000	−0.033078	0.000000	−0.033078	0.000000	−0.033078
6	0.010933	−0.007655	0.010933	−0.007655	0.010747	−0.007295	0.010933	−0.007655	0.010933	−0.007655	0.010933	−0.007655
7	0.010157	−0.008901	0.010157	−0.008901	0.010037	−0.008670	0.010157	−0.008901	0.010157	−0.008901	0.010157	−0.008901
8	0.015549	−0.012842	0.015549	−0.012842	0.015359	−0.012651	0.015549	−0.012842	0.015549	−0.012842	0.015549	−0.012842
9	0.024757	−0.019430	0.024757	−0.019430	0.024694	−0.019200	0.024757	−0.019430	0.024757	−0.019430	0.024757	−0.019430
10	0.040879	−0.031224	0.040879	−0.031224	0.041076	−0.031181	0.040879	−0.031224	0.040879	−0.031224	0.040879	−0.031224
11	0.020164	−0.005416	0.020164	−0.005416	0.019642	−0.005210	0.020164	−0.005416	0.020164	−0.005416	0.020164	−0.005416
12	0.020301	−0.004872	0.020301	−0.004872	0.019816	−0.004809	0.020301	−0.004872	0.020301	−0.004872	0.020301	−0.004872
13	0.030207	−0.008308	0.030207	−0.008308	0.029842	−0.008150	0.030207	−0.008308	0.030207	−0.008308	0.030207	−0.008308
14	0.046860	−0.013632	0.046860	−0.013632	0.046616	−0.013580	0.046860	−0.013632	0.046860	−0.013632	0.046860	−0.013632
15	0.100048	−0.028230	0.100048	−0.028230	0.100091	−0.028551	0.100048	−0.028230	0.100048	−0.028230	0.100048	−0.028230
16	0.028501	−0.002651	0.028501	−0.002651	0.028253	−0.002573	0.028501	−0.002651	0.028501	−0.002651	0.028501	−0.002651
17	0.029659	−0.001586	0.029659	−0.001586	0.029287	−0.001512	0.029659	−0.001586	0.029659	−0.001586	0.029659	−0.001586
18	0.041562	−0.003615	0.041562	−0.003615	0.041091	−0.003537	0.041562	−0.003615	0.041562	−0.003615	0.041562	−0.003615
19	0.062858	−0.006678	0.062858	−0.006678	0.062583	−0.006649	0.062858	−0.006678	0.062858	−0.006678	0.062858	−0.006678
20	0.103448	−0.010629	0.103448	−0.010629	0.103312	−0.010919	0.103448	−0.010629	0.103448	−0.010629	0.103448	−0.010629
21	0.030675	0.000000	0.030675	0.000000	0.030360	0.000000	0.030675	0.000000	0.030675	0.000000	0.030675	0.000000
22	0.033662	0.000000	0.033662	0.000000	0.033220	0.000000	0.033662	0.000000	0.033662	0.000000	0.033662	0.000000
23	0.046097	0.000000	0.046097	0.000000	0.045517	0.000000	0.046097	0.000000	0.046097	0.000000	0.046097	0.000000
24	0.068770	0.000000	0.068770	0.000000	0.068295	0.000000	0.068770	0.000000	0.068770	0.000000	0.068770	0.000000
25	0.105491	0.000000	0.105491	0.000000	0.105030	0.000000	0.105491	0.000000	0.105491	0.000000	0.105491	0.000000
Пp (Nmm)	−900.16		−900.16		−900.08		−900.16		−900.16		−900.16
Analysis number	59025		63115		59146		29997		56575		−

for FPA, JA, HS, HHS, and TLBO. Similar to Structure 1, HHS, FPA, JA, and TLBO have the same results with FEM, while HS is close to the solutions. The number of analyses to find the minimum energy is close to each other for all methods except HHS. In this example, one cannot make a distinction among the classical form algorithms applied to their accuracy, speed, and robustness, although the number of state variables becomes increased from 30 to 40 as compared to 1^st example. HSS also has a fast-computational ability comparing to all other algorithms.

4.3. Structure 3. A Plate with a Hole in the Middle (144 Members and 90 Nodes)

The plate presented as Structure 2 is also investigated with a different meshing option, as seen in Figure 7, in order to find a deflected form of the system for more points on the plane. Furthermore, by considering more members, an analysis result close to the exact solution can also be found. In that case, the total potential energy and the displacement values for this example will be different to the results of Structure 2. The minimum energy and the analysis number needed for the minimum value are given as

Table 3. Total potential energy results for the system with 144 members.

	TPOMA (FPA)	TPOMA (JA)	TPOMA (HS)	TPOMA (HHS)	TPOMA (TLBO)	FEM
Пp (Nmm)	−905.08	−905.08	−901.48	−905.08	−905.08	−905.08
Analysis number	999992	972853	876457	269416	947564	−

. All nodal displacements are plotted in Figure 8 for FEM and HHS results. Similar conclusions can be drawn in this example as compared to the first two problems solved.

5. Conclusions and Future Studies

The main purpose of this study is to see the applicability of the method called TPO/MA to more general structures than those it was previously applied to. In fact, TPO/MA has been applied until now to trusses and truss-like structures as cables and tensegric structures. The current study has shown that plates can be also analyzed using TPO/MA.

According to the No-Free-Lunch theorem [39], an algorithm may outperform another one in an example, while another one outperforms that algorithm in a second example. For that reason, several metaheuristic algorithms, such as HS, FPA, TLBO, and JA, were tested with the proposed novel hybrid method in the study. While HS has slightly failed to find the minimum energy level, the proposed HHS reached to minimum energy level quickly than the other algorithms. The methodology using HSS is effective on TPO/MA on solving plane stress problems. The processing time of the all iterations is in an acceptable value. By using TPO/MA, 1000 iterations of three examples were solved in 1.5 s, 2.5 s, and 8 s, respectively. Compared to other approaches, HHS can nearly find the global optimum results by half and 1/3 of the number of function evaluations for the first two examples and third example, respectively.

The analysis results found according to the HHS approach is effective to find the same results with the FEM method. Only minor differences (last than 1%) can be seen, but the total potential energy values are not affected.

Based on these considerations, it can be easily foreseen that this study will be followed by applications of TPO/MA on plane strain problems, plates with three-dimensional deformations, those with nonlinear constitutive equations, and those with other nonlinearities.