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

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.


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

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.
u(x, y) = u i + C 1 x + C 2 y, (1) v(x, y) = v i + C 3 x + C 4 y, 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 1 , C 2 , C 3 , and C 4 .
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: u(a j , b j ) = u j , v(a j , b j ) = v j , u(a k , b k ) = u k , v(a k , b k ) = v k , Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 16 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.

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 ui and vi 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.
When the particular derivatives of displacement fields are taken with respect to x and y, the strains are expressed in terms of the constants C1, C2, C3, and C4.
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:   (1) and (2), the relations between the nodal displacements is written as Equation (9): The solution of Equation (9) yields the constants C 1 -C 4 as in Equation (10)- (13): The strain energy density that resulted due to deformation is written as Equation (14) for a two-dimensional elastic body: When the solution of constants given as Equations (10)- (13) are used in Equations (3)-(5), the stains are written as follows: 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: The substitution of Equations (15)-(17) to (18)- (20) yields Equations (21)-(23) for stresses for linear problems: Then, the strain energy for linear problems can be found from Equation (14), as can be seen in Equation (24), for one element: 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): 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: 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 ):

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 ): 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. Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 16 In local optimization phase of HS, a randomly selected existing solution (Xran) is used, and the neighborhood range (Xran ± fw) is scanned. The formulation of the local optimization of HS is shown as Equation (31): The parameter "pitch adjusting rate (PAR)" is a parameter controlling the acceptable range around the chosen solution.

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

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.
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): The parameter "pitch adjusting rate (PAR)" is a parameter controlling the acceptable range around the chosen solution.

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.

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.

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 2 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.

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.

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 2 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, 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 cantilever beam has 15 free nodes and 30 displacement values considered as design variables. As seen in the results presented in Table 1, 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 6 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.

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 2 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 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.

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. 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. 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 2 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 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 2 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.    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. 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.

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

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.