Optimizing the buckling characteristics and weight of functionally graded circular plates using the multi-objective Pareto archived simulated annealing algorithm (PASA)

In this study for the first time, weight and critical buckling load in two kinds of functionally graded (FG) circular plates, namely, aluminum–alumina of (Al/Al2O3) and aluminum–zirconia (Al/ZnO2), are optimized using multi-objective Pareto archived simulated annealing algorithm (PASA). Material properties are assumed to vary with the power law in terms of the volume fractions of the constituent in two forms of symmetric and asymmetric with respect to the middle surface. The plate is subjected to uniform radial load and is considered for two boundary conditions, namely, simply supported and clamped edges. Aim at obtaining the Pareto archive is to achieve simultaneously the maximum buckling and the minimum weight concerning with proposed constraints. The parameters include the radius, thickness and volume fraction that the certain range is intended individually. The constraints are presented in form of the ratio of thickness to radius in category of the thin plates as well as the critical buckling stress being in the elastic range. Proposed simulated annealing algorithm is coded in MATLAB to obtain optimal non-dominated solution.


Introduction
Nowadays, the development of intelligent systems inspired the nature, is one of the most popular areas of artificial intelligence. Methods such as simulated annealing algorithms, genetic algorithms, ant colony optimization, are of the issues that human beings could achieve them by inspiring nature. On the other hand, plates are structures that are commonly used in military, petrochemical, aerospace industries, such as those found in turbine drive system, wall of pressure vessels, nuclear reactors and storage tanks floor. In recent years, functionally graded materials are one of the efficient and modern composite that are used in the abovementioned industry, especially for using in the high temperature environments.
These kinds of materials are advanced nonhomogeneous microstructure composites whose mechanical properties change smoothly and continuously from one surface to another. One of the most important phenomena in the plates, which causes the inefficiency in performance, is the buckling of the plate. In thin plates due to their low thickness, the membrane stiffness of the plate is much higher than its bending stiffness. Thus, it saves a large amount of membrane strain energy without causing much deformation. If the same amount of membrane energy is stored in the form of bending energy in the plate, it will be much deformed. Therefore, if the membrane strain energy stored in the plate is converted to bending strain energy under certain conditions, the plate is severely damaged that this phenomenon is called the buckling. Study of buckling behavior of plates have always been considered as one of the most important subjects in the structure analysis. Brush and Almorth [1] have comprehensively analyzed the buckling problems of columns, plates and shells and studied various methods for formulating the governing nonlinear equilibrium equations in buckling. In addition, many researchers have provided analytical and numerical solution for study of the buckling behavior of circular and annular FG plates based on thin, first and higher order shear deformation theories [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20].
In general, the optimization methods has been developed to improve the response to engineering demands in the modern industrial world; objective functions are selected depending to the performance criteria in engineering. For instance including weight loss, reduction of residual stresses, increase of vibration frequency, improvement of failure resistance, and increase of buckling load. So far, many studies have been carried out to optimize the behavior of the plates. Osaka et al. [21] examined the optimization of thickness profile of solid and annular plates in order to maximize critical buckling load using the finite element method and taking into account the shear deformation effects based on the first-order shear deformation. Ootao et al. [22] optimized the fraction volume index of FG plates by the neural network in the desired temperature range. In the other work, Ootao et al. [23] optimized the fraction volume index of FG hollow cylinder under thermal loading using neural networks method. Goupee and Vel [24] used the element-free Galerkin method and the genetic algorithm for two-dimensional optimization of the fraction volume index of functionally graded materials in two separate model problems. In one model, they minimized the residual stress by cooling down a FGM component (Ni/Al 2 O 3 ) and in another one, the weight of the Al/zirconia (FGM components) was optimized by restriction of the effective stress and maximum temperature experienced by metal constituent. Mozaffari et al. [25] optimized the critical buckling load of the FGM rectangular plate using the Imperialist Competitive Algorithm. Fereidoon et al. [26] used particle swarm-based algorithms to optimize the volume fraction of functionally graded materials. Juin-Der Duh and Brown [27] showed in multi-optimization problems, the solutions suggested by the knowledge-informed approach are more effective in approximation the set of Pareto optimal solution than those generated by the standard Pareto simulated annealing. Sahin and Turkbey [28] obtained Pareto solution for multi-objective facility layout using simulated annealing algorithm.
As abovementioned, buckling behavior of plates has been received the great deal of attention in scientific publications and is always considered as one of the most important subjects in structural analysis. Thus, the scholars have been paid attention to the optimization of buckling behavior. In the present study, two-objective optimization of buckling of the plates has been studied. The aim is to obtain the maximum critical buckling load according to the existing constraints with the approach of achieving the minimum weight. Thickening the plate to increase the critical buckling load, beside increasing the weight of the plate and consequently, more transverse loading, as well as increasing the cost of consumables materials. Therefore, it is important to determine the appropriate thickness of the plate and appropriate composition of volumetric fraction of ceramics and metal can achieve to the optimal bending rigidity required to maximize critical loads. Herein, the simulated annealing algorithm (SA) is exploited as the optimization method. The rationale behind SA extends to the annealing process of physical systems applied in thermodynamics. In this process, a physical system initially at a high-energy state is gradually cooled down until its minimum energy level is reached [29]. The idea is establishing a direct analogy between minimizing the energy level of a physical system and lowering the cost of an objective function. SA algorithm takes CPU time less than genetic algorithm (GA), due to find the optimal solution using point-by-point iteration rather than a search over a population of individuals [30]. The unique feature of SA accepts to move towards a worse solution during the search, therefore it prevents this method to trap in the local optima [31].
The aim of multi-objective optimization is to find a set of compromise solutions with different trade-offs among criteria, also known as Pareto optimal set; when this set is plotted in the objective space, it is called the Pareto front [32]. Serafini [33] presented the first version of the multiobjective simulation algorithm, was similar to the simulation algorithm of single objective optimization. The difference with the main algorithm was in the method of acceptance criteria. He investigated different alternative criteria to increase the probability of non-dominated optimal solutions that achieved to a special rule to obtain the non-dominated Pareto solutions. Suman [34] presented four simulated annealing based multi-objective algorithm and showed they are robust with algorithmic parameters and are capable of generating a large number of optimal solution. Suppapitnarm et al. [35] presented a novel implementation of the simulated annealing algorithm designed to explore the trade-off between multiple objectives. They suggested a new acceptance probability formulation based on proposition multiple temperature for multiple objectives in an annealing schedule. Czyzak and Jaszkiewicz [36] presented a new algorithm for multiobjective simulated annealing based on the Pareto optimal that is called PSA. Ulungu et al. [37] proposed and tested a multi-objective simulation algorithm (MOSA) to solve problems. They used a cumulative weight function to evaluate the objective function. The algorithm started with an initial solution, but a massive non-dominated solutions was obtained during the algorithm implementation. Akbulut and Fazil [38] optimized the thickness of a composite plate by simulated annealing algorithm (SA). They considered orientation of the resin and the number of sub-layers in each layer as the design variables and predicted failure. Moreover, Gutjahr and Pichler [39] used multi-dimensional optimization for studying non-computational methods as well as computational modeling to solve the problems. Recently, Rangaiah [40] described the process of engineering systems using multi-objective optimization procedure in his book. Chen et al. [41] compared the simulated annealing algorithm by using the tabu-search algorithms through study of advanced programming system for filtering color of TFT-LCD screens. Deb [42] provided a brief introduction and research focus to the field of multi-objective optimization based on evolutionary algorithm (EO). He remarked the multi-objective optimization methods give rise to set of Pareto-optimal solution, which by using the modified EO, the preferred region is provided on the Pareto-optimal front instead of complete front. In the engineering design community, size, shape, and topology optimization procedures are three classes of extensively studied design methodologies with assumed materials. Tang et al. [43] concerned with simultaneous designs of material selection and geometry optimization under static and thermal loads in a framework of multi-objective optimization of tracking the Pareto curve. Study of literature review reveals that this work is the first attempt to investigate the optimization of material composition and geometrical parameters of a FG circular plate by utilizing of a multi-objective simulation annealing algorithm using Pareto archive. The present paper is structured as follows. In Section 2, we give the expressions for the present problem, Section 3 is devoted to introduce Pareto archived simulated annealing algorithm. In Section 4, we present the numerical example for validation of the exploited algorithm. Finally, in Section 5, the performance of the proposed algorithm is investigated by studying the optimization in buckling and weight for two kinds of FG circular plate namely, aluminum-alumina of (Al/Al 2 O 3 ) and aluminum-zirconia (Al/ZnO 2 ).

Definition of the problem
A circular solid FG plate with radius R and thickness h under a radially uniform pressure loading is considered. The mechanical properties are assumed to be varied in the direction of the thickness. Figure 1 shows the schematic variation of the properties for symmetric and asymmetric FG plates that they can be defined as follows: where P(z), P m and P c denote the mechanical properties for the FG plate, the metal and the ceramic, respectively. f c denotes the volume fraction of the ceramic and varies between zero and one. Mechanical properties include elasticity modulus, density and yield stress of the FG material. Moreover, h and k are the thickness and the FG volume fraction index, respectively. As the critical buckling load and the weight of FG plates are depended to geometrical parameters and the volume fraction index of k, the purpose of this study is to determine the optimal values of these variables, so that the buckling load can be maximized, whereas the value of the weight is minimized. The relations for determination of the weight of solid circular plate and critical buckling load are as follows.
The weight of circular FG plate: The buckling load of circular plate [4]: (A) Critical buckling load for clamped supported edge: (B) Critical buckling load for simply-supported edge: In the above equations, R and D are the radius and the bending rigidity of the plate, respectively. For the asymmetric and symmetric FG plates, the bending rigidity is defined as follows: (A) Bending rigidity of symmetric the FG plate: (B) Bending rigidity of the asymmetric FG plate: In order to deal with problem of optimization of material and geometrical parameters to minimize weight and maximize critical buckling load, the Z-energy function is defined to estimate the metropolis criterion [44]: In the above relation, f 1 and f 2 stand for log(w) and log(P cr ), respectively. b denotes the weighting coefficient of the weight of plate, therefore (1Àb) evaluate the weighting coefficient of the critical buckling load. Noteworthy is that due to be not equal the scale of both objective functions, the logarithm of the two objective functions is utilized. Herein, b equal to 0.5, which means that both weight and buckling coefficient containing the identical significance for designer to make a decision among the optimal solutions according to the design task and its constraints. The aforementioned function is introduced as the sum of the two objective functions, and the critical buckling load is indicated by a negative sign in the energy function because the characteristics of the two objective functions is not the same, one function (plate weight) must be minimized and another one (critical buckling load) would be maximized.
The constraints are expressed as: where s Y is the mean yield stress. The latter constrain means that the design of the optimum plate is remained in the elastic region; to this end, critical buckling load should be less than s Y h. The problem is formulated as: By assuming that the ceramic constituent obeys linear elastic behavior, whereas the metal constituent achieves to the yield stress s Ym , with respect to the rule of volumetric mixture, the value of the mean yield stress at any arbitrary point is determined as follows [45]: It should be noted that the parameter q in equation (8) is determined numerically by the finite element micromechanics mode. Selecting a value of 500 GPa for q is appropriate in a wide range of volume fraction of FG plates with ceramic-metal composition [45]. Herein, the mechanical properties of metal-ceramic constituents of two kinds of FG plate, namely aluminum-zirconia (Al/ZnO 2 ) and aluminum-alumina (Al/Al 2 O 3 ) are given in Table 1.

Pareto archive of multi-objective annealing simulation algorithm
It is noteworthy to mention that contrary to the singleobjective optimization problem, for the multi-objective optimization problem, there is a set of non-dominated solution called Pareto optimal set. Simulated annealing (SA) is a stochastic search method, which initiates the physical annealing of solid for finding solution to combinatorial of annealing of solids [28]. This algorithm is one of the few algorithms that have explicit strategies to avoid local minima. For this reason, the solutions of worse quality than the current solution is accepted in order to escape from local minima [46]. The basic steps involved in the PASA algorithm for the present problem are as follows: (a) Define the simulated annealing parameters, such as: the initial temperature, cooling factor, weighting coefficients for two objective function. If it is met, the algorithm is finished. If the stopping criteria are not fulfilled, the algorithm goes to step b.
According with the aforementioned relationship, when the temperature is high, most of points are accepted so that the algorithm can search the entire target space and, with the progress of the algorithm, the probability of accepting the points that worsen the situation is reduced, and only the points that cause improvement of the algorithm are accepted. In the first step, after entering the initial parameters, initially the algorithm starts with a primary random point in the feasible range and evaluates this point, then the initial point is saved as the current point of the problem to Pareto archive as well as the optimal point. The parameters used in the simulation annealing algorithm are included of the initial temperature, the final temperature, the predefined number of iterations, the length of decrement, the relationship between the step of reducing the temperature and the initial solution. Herein, the initial temperature is determined by getting the average of 50 energy functions selected randomly T i ¼ Moreover, according to the stopping criteria of the algorithm, the final temperature is obtained using the equation of T f = 0.01T i . The number of iterations for the current temperature level chosen to be 200. Using the relation of the cooling function (T k+1 = aT k ), the temperature is decremented, in which a is the cooling coefficient, which is considered between zero and one. In this paper, a is chosen to be 0.98. The flowchart of the present method is described in Appendix A.

Providing a numerical example for validation
In this section, two objective functions are minimized multi-objective Pareto archived simulated annealing algorithm and compared the results with genetic algorithm using uniform genetic algorithm (UGA), hybrid genetic algorithm (HGA) [47][48]. The test objective functions are defined as follows [49]: So that -3 < x < 3, -5 < y < 5. In addition, Table 2 shows the number of generated solutions and the number of Pareto solutions using by different algorithms.  As shown in Table 2, the number of solutions obtained from the Pareto archive is more than the total evaluated solutions by the HGA algorithm. According to the results, the scattering of the results of the simulation annealing algorithm is greater, the optimal Pareto solution are searched in greater space.

Results
Herein, for the present problem, three approaches are represented. The results are demonstrated for the two metal-ceramic constituents of FG plate, as: aluminumzirconia (Al/ZnO 2 ) and aluminum-alumina (Al/Al 2 O 3 ).

The first approach: constant radius
In the first case, the radius is assumed to be constant, and the volume fraction and thickness of the plate is considered as variables that would be optimized. Figures 2-5 show the distribution of all dominated accepted and Pareto optimal solutions (Pareto optimal front) for asymmetric and symmetric FG plates with two boundary conditions, namely, simply-supported and clamped edges for radius of R = 1 m. In all figures, the bright blue and dark blue points indicate all the acceptable dominated and non-dominated optimal solutions (Pareto front), respectively. Figure 2 shows the results of the implementation of algorithm for Al/ZnO 2 plate with simply supported and clamped edges in the asymmetric and symmetric cases. When the radius of the plate is specified and the aim is to achieve optimum critical buckling load and weight, thus Figure 2 can be used. As sketched, firstly, the graph raises continuously, then a rough trend is occurred due to difference between the properties of two constituents of FG plate. In addition, in Figure 2, the results obtained by the present method are compared with those of Genetic algorithm. As seen, the results indicate a very good agreement between them.
As illustrated in Figure 3, at the left side of boundary of the thickness, the optimal dark points indicate to the maximum volume fraction and the minimum value of the weight of the plate. In Figure 3a, while being far from the left side border, Pareto archive converge to a volume fraction index of 10-50 (that means a FG plate with a higher metal percentage). This trend continues until toward the end of the border on the horizontal axis labeled thickness. These points are not dominated by the greater value of thickness and the volume fraction index, due to the fact that the thickness limitation does not permit to proceed more; hence the volume fraction index tends to be zero to maximize the buckling load and minimize the weight. Those points are not added into Pareto archive if they violate the elasticity relation for calculation critical buckling forces.
The results of the algorithm are presented for Al/ Al 2 O 3 plate with a simply supported and clamped edges in asymmetric and symmetric cases in Figures 4 and 5.  mechanical properties of aluminum and alumina, as well as aluminum and zirconia, there is no convergence of the algorithm towards pure metal, and the optimal Pareto points are in the form of ceramic and metal mixture and tend towards ceramic properties. As sketched in Figure 5, at the initiation and the end of the thickness boundary, the volume fraction index converges to the upper boundary of the volume fraction index (pure metal). This is due to the fact that the weight of the plate in the boundary thickness is minimized at the beginning and end of the thickness boundary with high volume fraction index; these points do not violate clause i of the algorithm, therefore are added to the Pareto archive and is not dominated by any point.
In the upper part of Figure 4c and d, some dispersed distribution of Pareto optimal points is seen, this is due to achieving to the upper boundary of the thickness. As the constraint for thickness does not permit to be increased this parameter in the upper bound, the volume fraction index increases and if it does not violate part I of the algorithm, it is added to the Pareto archive due to containing the lowest weight. This can be observed for the corresponding points in diagrams (c) and (d) of Figure 5.

Second approach: constant volume fraction index
In this case, the volume fraction index is assumed to be constant, whereas the radius and thickness are variable. The range of variables in this case is 0.0125 < h R < 0.1 and 0.1 < R < 1. In this approach, the purpose is to estimate an optimum ratio of thickness to radius in the desired volume fraction index. Figures 6 and 7 show the results of an optimal ratio of thickness to radius for (Al/Al 2 O 3 ) and (Al/ZnO 2 ) FG plates corresponding to the prescribed volume fraction index, respectively. Through Figures 6 and 7, the volume fraction index is assumed to be taken of 0, 0.5, 2 and 10 for the optimal ratio of thickness to radius. For instance, if the volume fraction index of 10 is selected in Figure 7d, thus h R would be 0.061824 to be achieved to the maximum critical buckling load and the lowest weight. In the following, Figure 8 is sketched to present the Pareto archive corresponding to Figure 7d for aforementioned quantities of volume fraction index and thickness to radius ratio.
In Figure 8, the brighter area represents all the acceptable points and the darker curved lines represent the Pareto archive of non-dominated solution, which are   Figures 7 and 8, as the volume fraction index increases, the critical buckling load decreases; whereas the optimal ratio of thickness to radius increases; which consequently raises the weight of the plate.

Third approach: variable radius, thickness and volume fraction
In this case, the radius, thickness, and volume fraction of metal-ceramic constituents, are assumed to be varied. For instance, as presented in Tables 3 and 4, the weight of 100 kg is chosen as the target and the closest nondominated points for aluminum-zirconia and aluminumalumina plates are presented in Tables 2 and 3, respectively.
By running computer code in Matlab for three times, the solutions are tabulated in Tables 3 and 4. As shown, the objective is to obtain the closest non-dominated point to the specified weight; in fact, the obtained points containing the maximum buckling load. On the other hand, the aim can be considered as the closest non-dominated solution to a specified buckling load. In Tables 5-8, according to the buckling load of 5 Â 10 6 (N/m), the nearest non-dominated Pareto point is represented.  In the present article, the critical buckling load as well as the weight of the plate was optimized by the multiobjective annealing simulation algorithm using the Pareto archive. The aim is to find a non-dominated Pareto archive, so that the critical buckling load is maximized whereas the amount of plate weight is minimized. The optimization of the asymmetric and symmetric functionally graded aluminum-zirconia (Al/ZnO 2 ) and aluminum-alumina (Al/Al 2 O 3 ) plates under simply supported and clamped edge conditions have been considered. As the mechanical properties of traditional functionally graded materials depend on the volume fraction index of ceramic-metal composition, it is noteworthy to mention that which volume fraction index and thickness can satisfy all the constraints of the problem, and optimize the critical buckling load and weight. Herein, the solution based on the multi-objective annealing simulation algorithm is obtained using programming in MATLAB. The results are proceeded in three approaches and the solutions are extracted in three categories, as exhibited in tables and figures. The solutions are demonstrated as the non-dominated Pareto archive as well as the optimal solutions of the problem, In fact, it provides the decision-maker with a Pareto solution set and permit her or him to make a choice depending on design preferences. The closest and best solution of the Pareto archive solutions can be selected by accessing the Pareto archive and having enough information about it. Herein, the following conclusion can be outlined: I. By assuming that the maximum weight/ minimum critical buckling load is specified, the second approach can be used to obtain optimal state. II. By considering the necessities of the problem and using the third approach, and with the implementation of the algorithm for several times, the solutions based on the other two approaches can be obtained. III. As shown for Al/ZnO 2 plate, optimization of critical buckling load and weight is achieved from a plate containing a higher volume percentage of metal constituent; whereas, for Al/Al 2 O 3 plate, this trend is reverted and is resulted from a plate containing a higher volume percentage of ceramic constituent. IV. The present result is utilized as a new reference for the designers in the field of aerospace and marine structures industries. For future studies, optimization of critical temperature of thermal buckling, critical frequency of vibration of FG plate are interesting subjects to be considered.