Multi-Objective Optimization of Functionally Graded Beams Using a Genetic Algorithm with Non-Dominated Sorting

A mixed layer-wise (LW) higher-order shear deformation theory (HSDT) is developed for the thermal buckling analysis of simply-supported, functionally graded (FG) beams subjected to a uniform temperature change. The material properties of the FG beam are assumed to be dependent on the thickness and temperature variables, and the effective material properties are estimated using either the rule of mixtures or the Mori–Tanaka scheme. The results shown in the numerical examples indicate the mixed LW HSDT solutions for critical temperature change parameters are in excellent agreement with the accurate solutions available in the literature. A multi-objective optimization of FG beams is presented to maximize the critical temperature change parameters and to minimize their total mass using a non-dominated sorting-based genetic algorithm. Some specific forms for the volume fractions of the constituents of the FG beam are assumed in advance, such as the oneand three-parameter power-law functions. The former is used in the thermal buckling analysis of the FG beams for comparison purposes, and the latter is used in their optimal design.


Introduction
Functionally graded (FG) beams, plates, and shells are emerging heterogeneous material structures, which are formed by mixing two-or multiple-phase materials with a pre-designed spatial distribution of volume fractions of the constituents [1][2][3][4]. Because the material properties of FG structures gradually and smoothly vary through their domain, some drawbacks can be prevented, including delamination and stress concentration, which usually occur at the interfaces between adjacent layers in the case of the laminated composite structures due to the material properties suddenly changing at these locations. On the other hand, FG structures provide a large optimal design space in engineering practice. Engineers can design the spatial distributions of the volume fractions of the constituents according to practice demands to obtain the best structural physical properties. Based on the above-mentioned advantages of FG structures, they are also becoming increasingly more popular in various high-end industries, such as aerospace, submarine, automobile, and nuclear industries, as well as for use in biomedical implants [5,6]. A variety of mechanical analyses of FG structures and single/multi-objective optimization of these structures, thus, are attracting considerable attention.
Some comprehensive literature surveys with regard to articles examining the mechanical analyses of FG structures can be found in the public literature [7][8][9][10][11][12][13][14]. Among these, the literature survey in this work focuses on articles investigating the thermal buckling behavior of FG beams using various beam theories and on articles related to the singleand multi-objective optimizations of FG beams.
Based on the Euler-Bernoulli theory (EBT), Kiani and Eslami [15] examined the thermal buckling behavior of FG beams under different boundary conditions, in which only the bending effects were considered in their analysis, rather than the shear deformation effects. The first-order shear deformation theory (FSDT) was used by Tauchert [16] to study the thermal buckling of simply supported, antisymmetric angle-ply laminates, in which the shear deformations of the laminates were simply assumed as a constant through their thickness direction. Reddy [17] developed a refined shear deformation theory (RSDT) to investigate the various mechanical behaviors of laminated composite (LC) structures. In Reddy's formulation, the shear deformations induced in the LC structures were assumed to be a parabolic function distribution through the thickness direction; the shear correction factor commonly used in the FSDT was not needed, and the surface conditions regarding the traction stresses free were satisfied. In order to more precisely capture the shear deformation effects on the various mechanical behaviors of LC and FG structures, some advanced shear deformation theories have been proposed, including the third-order shear deformation theory (TSDT) [18,19], the sinusoidal shear deformation theory (SSDT) [20,21], the hyperbolic refined shear deformation theory (HRSDT) [22,23], and the exponential shear deformation theory (ESDT) [24]. Liu et al. [25] investigated the nonlinear response of FG structures using an iso-geometric continuum shell element method. Liu et al. [26] presented an analysis of thermo-mechanical snap-through instability of structures using a hybrid load-and displacement-controlled algorithm. Based on a shear deformation theory, Refrafi et al. [27] examined the hygro-thermo-mechanical buckling behavior of simply-supported FG sandwich plates resting on a Winkler-Pasternak foundation. Based on Carrera's unified formulation (CUF) [28], the above-mentioned advanced and refined shear deformation theories can be hierarchically derived and included as special cases [29]. Some comparative studies for different mechanical behavior analyses of FG beams using assorted advanced and refined shear deformation beam theories have been carried out, including thermal buckling and post-buckling analyses [30], bending and free vibration analyses [31], thermal buckling analysis [32], and free vibration analyses [33,34]. The above-mentioned theories can be summarized as equivalent-single-layered (ESL) theories, in which a global coordinate system was introduced to describe the displacement field induced in LC and FG structures.
Because ESL theories are often inadequate for determining the three-dimensional (3D) stress field induced in loaded LC and FG beams, some layer-wise (LW) theories for these beams have been presented. Based on Hamilton's theory, Tahani [35] developed a displacement-based LW theory for the static bending and free vibration analyses of LC beams. Lee and Saravanos [36] developed an LW FSDT for coupled thermo-piezoelectric composite beams under thermal loads. Shimpi and Ainapure [37] and Shimpi and Ghugal [38] developed an LW trigonometric shear deformation theory for the static analysis of LC beams. Based on an LW higher-order shear deformation theory (HSDT), Pandey and Pradyumna [39] developed a finite element method (FEM) for the static and dynamic analyses of FG sandwich plates. Bayat and Toussi [40] carried out thermal buckling and thermal postbuckling analyses of LC beams reinforced with shape memory alloy wires using a displacement-based LW theory. Based on the stationary principle of minimum potential energy [41,42] combined with the Lagrange multiplier method, Wu and Kuo [43] and Wu and Chen [44] developed a mixed LW HSDT for the bending, free vibration, and buckling behavior of LC plates. Wu and Xu [45] extended this approach to develop the strong and weak formulations of a mixed higher-order shear deformation theory for the static analysis of FG beams under thermo-mechanical loads. A comprehensive overview with regard to the development, numerical implementation, and application of LW theories was undertaken by Liew et al. [46]. Some articles related to the optimal design of LC structures and FG beams have been presented using different single-and multi-objective optimization algorithms combined with advanced and refined shear deformation beam theories. Walker et al. [47,48] developed an optimal procedure to select the best material combinations for minimum weight and cost of hybrid composite plates and for minimum mass of sandwiched composite cylindrical shells, respectively. Houmat [49] and Guenanou and Houmat [50] proposed a layer-wise optimization method to obtain an optimal lay-up design for maximum fun-damental frequency of variable stiffness LC plates and symmetrically LC circular plates, respectively. Cho and Ha [51] presented volume fraction optimization for the purpose of minimizing thermal stresses induced in Ni-Al 2 O 3 FG beams. Based on the Vlasov thin-walled theory [52], Nguyen and Lee [53] committed to the optimal design of thinwalled FG beams for buckling problems, in which an LW cubic interpolation for the through-thickness distribution of the volume fraction of the constituents was assumed, and a genetic algorithm (GA) was used as the optimal tool [54,55]. Kamarian et al. [56] studied the optimization of material compositions of FG beams to maximize their fundamental natural frequencies using the firefly algorithm and an adaptive neuro-fuzzy inference system. Yas et al. [57,58] engaged in the optimization of material compositions of a three-parameter FG beam in order to maximize its fundamental frequencies using the imperialist competitive algorithm, an artificial neural network [59], and a generalized differential quadrature (DQ) method [60][61][62]. Goupee and Vel [63] proposed a methodology to maximize the fundamental frequencies of FG structures by tailoring their material distributions. Na and Kim [64][65][66] investigated the volume fraction optimization of FG plates and panels under thermal loads to minimize the induced thermal stresses and to maximize the critical temperature parameters, in which a bi-objective optimization was considered and the optimal tool was a quasi-Newton method. Walker and Smith [67] presented a multi-objective optimization technique for minimizing the weighted sum of the mass and deflection of LC structures, in which a GA and an FEM were used. Tornabene and Ceruti [68] dealt with mixed static and dynamic optimization of four-parameter FG structures. The FSDT combined with the generalized DQ method was used for the static and dynamic analyses of the FG structures, and three different optimization schemes were used as the optimal tools, including the particle swarm optimization, the Monte Carlo, and the GA approaches.
In this work, the authors aim at investigating the material composition optimization of an FG beam under a uniform temperature change in order to maximize the critical temperature change parameter and minimize the total mass of the FG beam. A mixed LW HSDT [45] is developed for the thermal buckling analysis of the FG beam. The current solutions for the critical temperature change parameters of the FG beam obtained through considering the temperature-dependent (TD) and the temperature-independent (TI) material properties are examined, and solutions obtained using the von Kármán geometric nonlinear strains (GNS) and the full GNS are also estimated. A non-dominated sorting-based GA [69,70] is used for the current multi-objective optimization analysis in the current issue. The through-thickness distribution of the material properties of the FG beam is assumed to be a three-parameter power-law function of the volume fractions of the constituents [71,72], the material-property gradient indices of which are, thus, to be determined for the optimal material profile.

Effective Material Properties
In this work, two micromechanics models, the rule of mixtures [73,74] and the Mori-Tanaka scheme [75], are used to estimate the effective material properties of the FG beam, as described in the following sections.

The Rule of Mixtures
According to the rule of mixtures, the through-thickness distributions of the effective material properties of the FG beam are written in the following form, where Γ c and Γ m represent the volume fractions of the ceramic and metal materials of the constituents of the FG beam, respectively, such that Γ c + Γ m = 1. F can be one of the engineering constants, including Young's modulus E, Poisson's ratio υ, thermal expansion coefficient α, and mass density ρ. The subscripts m and c represent the metal material and the ceramic material, respectively.

The Mori-Tanaka Scheme
Within the frame work of Eshelby's elastic inclusion theory [73,74] restricted to a single inclusion in a semi-infinite elastic, homogeneous and isotropic medium, Mori and Tanaka [75] extended the original theory to more general cases of multiple inclusions embedded into a finite elastic medium, in which an equivalent inclusion-average stress method was used.
For a two-phase isotropic material, the effective material properties are written as the explicit formula [76,77] as follows: , and K and G represent the volume and the shear moduli, respectively. In later sections of this paper, if no specific mention is given the effective material properties of the FG beam are estimated using the rule of mixtures.

The Mixed LW HSDT for FG Beams
In this work, the authors develop a unified formulation of the mixed LW HSDT for the thermal buckling analysis of simply supported FG beams subjected to a uniform temperature change, in which the material properties of the FG beam are considered to be temperature-dependent. The configuration and coordinates of the FG beam are shown in Figure 1a, in which h and L represent the thickness and the length of the FG beam, respectively. In the analysis, the FG beam is artificially divided into N l layers as shown in Figure 1b

Strong-Form Formulation
In the current LW HSDT, the displacement field for a typical individual layer is given as follows: where m = 1, 2, . . . , N l ; u According to the perfect bonding assumptions at the interfaces between adjacent layers, the corresponding displacement continuity conditions at these places are given as where k = 1, 2, . . . , (N l − 1). The strain-displacement relationship is given as where the commas denote the derivative of the suffix variable, and remaining strains are zeroes, including ε . It is noted that the displacement continuity conditions at the interfaces between adjacent layers are imposed in this formulation, while the strain continuity conditions at these places are discontinuous. The stress-strain relationship for an orthotropic material in thermal environment is given as where the symbol ∆T denotes the temperature change, measured from a room temper- 3 , in which the subscripts 1, 2, and 3 denote the principle axes of the material properties, and that E, υ, and α represent the Young's modulus, Poisson's ratio, and the thermal expansion coefficients. For isotropic materials, these stiffness coefficients will be reduced to Q (m) Note that these engineering constants E, υ, and α in the analysis are considered to be dependent from the thickness and the temperature variables, and ∆T is a constant. In addition, other stress components are xy , because a two-dimensional beam model is used. The governing equations and associated boundary conditions are derived using the stationary principle of minimum potential energy combined with the Lagrange multiplier method, in which the displacement continuity conditions at the interfaces between adjacent layers given in Equations (7) and (8) are multiplied by the Lagrange multipliers and then substituted into the potential energy functional as the constraints, such that the extended potential energy functional of the N l -layered FG beam is given, as follows: where λ (m) x and λ Applying the stationary principle of minimum potential energy, following a standard variational process, and integrating the stress variables through the thickness direction, the authors obtain where Performing Equation (14) the integration by part yields According to Equation (15), the Euler-Lagrange equations of the mixed LW HSDT can be obtained as follows: δ w (m) δ w (m) δ w (m) where m = 1, 2, . . . , N l and k = 1, 2, . . . , (N l − 1). The possible boundary conditions are given as Either P Either R (m) Either N (m) Either M Either P (m) Either R where u Substituting Equations (A1)-(A11) into Equations (16)- (25), the authors obtain the Euler-Lagrange equations of the current mixed LW HSDBT in terms of all displacement components, and they are given as follows: δ w (m) δ w (m) where m = 1, 2, . . . , N l and k = 1, 2, . . . , (N l − 1). The total number of Euler-Lagrange equations, i.e., Equations (34)- (43), is (10N l − 2), (8N l − 2), and (6N l − 2) for the mixed layer-wise third-order (LW3), second-order (LW2), and first-order (LW1) shear deformation theories, respectively, are in terms of the same num-ber of unknowns as those of the corresponding Euler-Lagrange equations. The strong-form formulations of the mixed LW HSDT are, thus, obtained, including the Euler-Lagrange equations (Equations (34)- (43)) and the possible boundary conditions (Equations (26)-(33)).

Applications
Equations (34)-(43) associated with a set of boundary conditions (i.e., Equations (26)-(33)) can be composed as a well-post boundary value problem, and its Navier-type analytical solutions for the thermal buckling behavior of simply supported, multi-layered FG beams can be obtained using the Fourier series expansion method.

Statement of the Optimization Problem
In this work, the authors consider the multi-objective optimization of the volume fractions of the constituents of FG beams subjected to a uniform temperature change in order to maximize the critical temperature change parameter of the FG beam and to minimize its total mass. The FG beam is considered to be formed by mixing a metal material and a ceramic material according to a three-parameter power-law function of the volume fractions of the constituents through the thickness coordinate of the FG beam, which is given as follows: The three-parameter power-law function, where the symbols κ b , κ c , and κ p denote the material-property gradient indices. The material at the top surface of the FG beam (i.e., z = h/2) is ceramic rich, and when κ b = 0, the material at the bottom surface of the FG beam (i.e., z = −h/2) is metal rich. Variations in the through-thickness distributions of the volume fraction given in Equation (58) with varying one of κ b , κ c , and κ p and holding the other two indices constants are shown in Figure 2. In addition, when κ b = 0, the three-parameter power-law function is reduced to the one-parameter power-law function, which is commonly used in the literature and is given as follows: One-parameter power-law function, which is used in the thermal buckling cases in this work for a comparative study. Variations in the through-thickness distributions of the volume fraction given in Equation (59) with varying the κ p value are shown in Figure 3.  In the optimal design, a mass ratio R m is defined as follows: whereρ f ,ρ c , andρ m are the total mass per unit area in x-y plane of the FG beam, the homogeneous ceramic material, and the homogeneous metal material, respectively.
in which k = f, c, or m. The critical temperature parameter ∆T cr is defined as follows: where α c denotes the thermal expansion coefficient of a reference ceramic material.
The objective functions are defined as follows: Objective function 1 : Objective function 2 : where the ranges of F 1 and F 2 are 0 ≤ F 1 (or F 2 ) ≤ 1.
Because non-dominated sorting is used in the current GA, the values of these objective functions given in Equations (62) and (63) are, thus, used to classify each design into its corresponding non-dominated front, which is also the assigned fitness value used for the sorting process. Minimization of the fitness function can be accomplished when the critical temperature change parameters of the FG beam are obtained. The detailed process of the current non-dominated sorting-based GA is described in the following sections.

The Non-Dominated Sorting-Based GA
In this work, the non-dominated sorting-based GA [69,70] is used as an optimal technique, the usual form of which is described as follows: The GA starts with an initial set of random solutions, which are called a population. Each individual in the population is called a chromosome, which represents a solution to the problem at hand. The chromosomes are made of discrete units called genes (or design variables). Each gene controls one or more features of the chromosome. The chromosomes evolve through successive iterations, called generations. Based on some measures of fitness, the chromosomes during each generation are evaluated. To create the next generation, new chromosomes in the next generation, called offspring, are formed using two operators, which are the crossover and mutation operators. In the former, some portions of two chromosomes in the current generation are merged together, and in the latter, some portions of chromosomes in the current generation are modified. Thus, the crossover operator leads the population to converge by means of making the chromosome in the population alike, and the mutation operator assists the search escaping from local optima by reintroducing genetic diversity back into the population. A new generation is formed by selecting and rejecting some of the parents and offspring according to their fitness values so as to keep the population size of each generation constant. After several generations, the GA converges to the best chromosome, representing the optimum solution to the problem considered. A flow chart of the non-dominated sorting-based GA is shown in Figure 4, and its related process is described as follows: (a) Initial population The first population consists of chromosomes (i.e., solutions or designs), which is taken to be 200 in this work. Each chromosome is formed by several genes (i.e., design variables), each of which is represented by a randomly generated real number with four decimal places in a specific range.

(b) Fitness value
As mentioned before, there are two objective functions considered in this work, the critical temperature change parameter of the FG beam and its total mass, the specific forms of which (i.e., F1 and F2) are, thus, defined as the fitness values and are used for the sorting process in the current GA.

(c) Selection
In this work, the selection process consists of two approaches, the Pareto-ranking approach [55] and the crowding distance approach [69], which are described as follows: (c-1) Pareto-ranking approach The Pareto-ranking approach uses the Pareto dominance concept to evaluate fitness or assign selection probability to solutions. Each solution in the population of the parents is ranked according to a dominance rule, and then each solution is assigned a fitness value based on its rank in the population, rather than its actual objective function value. A dominance rule proposed by Goldberg [55] is adopted and is described as follows: Step 1: Set i = 1 and n p 0 = the number of original solutions in the population.
Step 2: Identify non-dominated solutions x j (j = 1~n nd f i ) in n p i , and assign them set to NDF i , in which NDF i denotes the ith non-dominated front, and the number of solutions in NDF i is n nd f i .
Step 3: Set i = i + 1, and then n p i+1 = n p i−1 − n p i . If n p i+1 = 0, go to Step 4; otherwise, go to Step 2.
Step 4: For every solution x k , k = 1 ∼ n p 0 at generation t, assign rank In the above process, a solution A dominates the other solution B, which means both Because all objective functions are to be minimized, a lower rank corresponds to a better solution, such that the population will be classified as several non-dominated fronts, where the first non-dominated front is also called the Pareto fronts of this population.
(c-2) Crowding distance approach A crowding distance approach as proposed by Deb [69] is adopted in this work, which is aimed toward obtaining a uniform spread of solutions along the best-known Pareto front. The approach is described as follows: Step 1: Rank the population and identify Pareto fronts as PF 1 , PF 2 , . . . , PF R , using Equations (62) and (63). For each non-dominated front, repeat Steps 2 and 3.
Step 2: For each objective function F k , sort the solutions in PF j in the ascending order. If the number of solutions in PF j is l, then the crowding distance of the ith solution is measured as where i = 2, 3, . . . , (l−1), and cd k (x 1 ) = cd k (x l ) = ∞.
Step 3: To obtain the total crowding distance cd( The selection rule is based on the rank assigned to the solutions, where a lower rank corresponds to a better solution. When the solutions are in the same rank, a greater crowding distance corresponds to a better solution.

(d) Crossover
Any two solutions in the population are randomly paired for mating, in which a simulated binary crossover operator suggested by Deb and Goyal [78] is adopted in this work and is given as follows: x j t+1 2 where γ j = 2u j 1/(η j +1) when u j ≤ 0.5, and γ j = 1/ 2 1 − u j 1/(η j +1) when u j ≥ 0.5, in which u i denotes a randomly generated real number between 0 and 1 and η > 0 is a user-defined index parameter, the value of which is suggested to be η j = 1.

(e) Mutation
Mutation occurs in each solution in the population when the randomly generated number between 0 and 1 for each solution is less than 0.1. A polynomial mutation operator suggested by Deb and Goyal [78] is adopted in this work and is given as follows: where in which u i denotes a randomly generated real number between 0 and 1 and η i is a userdefined index parameter, for which the value is suggested to be η i = 20.
(f) Elitist Performing the above operators, including the non-dominated sorting, crowding distance, crossover, and mutation approaches, we can obtain the new generation (i.e., offspring) with the same number of solutions (i.e., n p ) as those in the original generation (i.e., parents). The elitist strategy in this work is as follows: Step 1: Combine the solutions of parents and offspring to form a population with 2 n p solutions.
Step 2: Perform the non-dominated sorting and crowding distance approaches to each solution of the population generated in Step 1. Then, select the first n p solutions with smaller fitness values in the population to form the offspring.

(g) End criteria
The end criteria of the current non-dominated sorting-based GA is given as Criteria: when the number of generations reaches an assigned number, which is 100 in this work.

Thermal Buckling Analysis of Laminated Composite Beams and FG Beams
Because no benchmark solutions for the critical temperature changes for a simplysupported FG beam can be found in the literature, for comparison purposes, the current mixed LW HSDT is applied to the thermal buckling analysis of simply-supported, cross-ply laminated beams under a uniform temperature, in which the material properties are regarded as a layer-wise constant distribution varying through the thickness coordinate. The current solutions for critical temperature change parameters are compared with those obtained using the EBT [79], FSDT [79], TSDT [79], HSDT [80], and HRSDT [22]. The material properties and geometric parameters are given as L/h= 10 and 50, E 1 /E 2 = 10, 20, and 30, G 12 /E 2 = G 13 /E 2 = 0.6, υ 12 = υ 13 = 0.25, α 2 /α 1 = 3, in which the subscript 1 represents the material properties parallel to the reinforced fiber direction, and the subscripts 2 and 3 represent material properties perpendicular to the reinforced fiber direction. A dimensionless temperature change parameter is defined as ∆T cr = (∆T cr ) α 1 (L/h) 2 . Tables 1 and 2 show the convergent study for critical temperature change parameter solutions of [0 • /90 • /0 • ] laminated beams subjected to a uniform temperature change. It can be seen in Tables 1 and 2 that the current LW1, LW2, and LW3 solutions converge rapidly. For a moderately thick beam (L/h = 10), the convergent solutions are yielded when the total number of layers are taken to be 9, 6, and 3 for the LW1, LW2, and LW3 theories, respectively, and these convergent solutions closely agree with one another. The convergence rate increases when the beam becomes thinner. The current solutions in Tables 1 and 2 with a superscript a indicate that the von Kármán GNS, rather than the full GNS, is considered. The results show that the deviations between the solutions considering the von Kármán GNS and the full GNS are less than 1%, which is very minor, such that the von Kármán GNS is suitable for the current analysis of moderately thick laminated composite beams. The current solutions in Tables 1 and 2 with a superscript b indicate that a narrow beam is considered, in which Poisson's ratios are taken to be zeroes in the stiffness coefficients Q ij (i, j = 1-3). The results show that the deviations between the solutions considering and those not-considering the narrow beam effects are significant, where the deviations as high as approximately to 5.5% for both a moderately thick (L/h = 10) and a thin LC beam (L/h = 50). Tables 1 and 2 also show that the current convergent solutions considering the narrow beam effects closely agree with the HSRDT [22] and TSDT [79] solutions, and those not considering the narrow beam effects closely agree with the HSDT solutions [80]. In addition, the EBT solutions for the cases of L/h = 50 are the same as those for the cases of L/h = 10. In the former, the relative error between the EBT solutions and the current LW HSDT solutions is 1.2%, while it is 31.2% for the latter, such that the EBT is suitable for the current analysis of very thin laminated composite beams only. Table 1. Convergent study of critical temperature parameter solutions of simply supported, laminated (0 • /90 • /0 • ) beams under a uniform temperature change (α 2 /α 1 = 3, E 1 /E 2 = 20).
The thermal buckling problems of simply-supported FG beams with either the TD or the TI material properties under a uniform temperature change are considered in Table 3. A flow chart for the solution process to determine the critical temperature change parameters of the FG beam is given in Figure 5 where the TD material properties are considered, in which the golden section method is used. The through-thickness distributions of the material properties of the FG beam are assumed to obey a single-parameter power-law function according to the volume fractions of the constituents, as shown in Equation (57).
The effective material properties are estimated using either the rule of mixtures or the Mori-Tanaka scheme, which are given in Equations (1) and (2), respectively. Again, it can be seen in Table 3 that the current LW HSDT solutions converge rapidly. The solutions for the critical temperature change parameters obtained using the Mori-Tanaka scheme are approximately no more than 2% greater than those obtained using the rule of mixtures. The effects of different micromechanics models on the critical temperature change parameters of the FG beam are minor. Table 2. Convergent study of critical temperature parameter solutions of simply supported, laminated (0 • /90 • /0 • ) beams under a uniform temperature change (α 2 /α 1 = 3, L/h= 10). 0.8965 HSDT [80] 0.81683 HRSDT [22] 0.8833 TSDT [79] 0.8832 FSDT [79] 0.8868 EBT [79] 1. 0.7720 HSDT [80] 0.72608 HRSDT [22] 0.7472 TSDT [79] 0.7471 FSDT [79] 0.7528 EBT [79] 1.1329 a von Karman geometric nonlinear terms are considered. b The Poisson ratio is considered to be zero in Q ij (i = 1 and 3).
The results also show that the critical temperature change parameters obtained using the TI material properties always are greater than those obtained using the TD material properties. For a moderately thick beam (L/h = 10 and κ p = 1), deviations between the critical temperature change parameters obtained using the TI and the TD material properties are as high as approximately 36.5%, such that the current analysis using TI material properties is an unsafe analysis and, thus, is not recommended. The effects of TD material properties on the critical temperature change parameters are significant, and these must be considered for the current issue.

Optimization of Material Composition of FG beams
In this section, a multiple-objective optimization for material composition of a simplysupported FG beam subjected to a uniform temperature change is considered in order to maximize its critical temperature change and minimize its self-weight. In the optimal design, the FG beam is considered to be a two-phase composite material, where one phase is the ceramic material (ZrO 2 ) and the other is the metal material (SUS304). The TD material properties of ZrO 2 and SUS304 are given in Table 4 [4]. Variations in the material properties of ZrO 2 and SUS304 with the temperature variable ranging from 300 K to 1100 K are shown in Figure 6. The material properties of the FG beam are assumed to obey a three-parameter power-law distribution of volume fractions of the constituents along the thickness of the FG beam, and the effective material properties are estimated using the rule of mixtures. The non-dominated sorting-based GA mentioned in Section 4.2 is used to determine some sets of Pareto-optimal solutions of the undetermined coefficients κ p , κ b , and κ c , in which 200 initial populations are randomly generated, in which the ranges of κ p , κ b , and κ c are taken as0 ≤ κ p ≤ 50, 0 ≤ κ b ≤ 1, and 1 ≤ κ c ≤ 3, respectively.   8 show the populations at the initial, first, second, fifth, tenth, and 20th generations, where the TI and TD material properties are considered, respectively. It can be seen in Figures 7 and 8 that the non-dominated, sorting-based GA converges rapidly and that the Poreto-optimal solutions can be yielded after 20 generations. Curve fitting for these Preto-optimal solutions is performed for the cases of TI and TD material properties, which is shown in Figure 9 with the solid and dashed lines, respectively, in which 9936for the TD material property cases.
Equations (68) and (69) can then be used to calculate the associated weight numbers (i.e., w 1 and w 2 ) for the objective functions F 1 and F 2 , in which if we let dF 2 /dF 1 = λ, then w 1 = −λ/(1 − λ) and w 2 = 1/(1 − λ). Tables 5 and 6 show every ten Pareto-optimal solutions sorted by the F 2 function values from the largest to the smallest for the TI and TD material property cases, respectively. These Pareto-optimal solutions may provide design engineers with valuable information regarding what set of the material-property gradient indices (κ p , κ b , and κ c ) they need according to the weight number ratio of (w 2 /w 1 ). In addition, the results also show that the critical temperature changes of the FG beam for the TD material property cases are much less than those for the TI material property cases, such that the TD material property effects must be considered in TD problems.   Curve fitting for the Pareto-optimal fronts for the cases of TI and TD material properties. Table 5. Pareto-optimal solutions for simply supported FG beams under a uniform temperature change sorted according to the F 2 function values from the largest to the smallest, in which temperature-independent material properties are considered.

Conclusions
In this work, the authors developed a mixed LW HSDT for the thermal buckling analysis of FG beams subjected to a uniform temperature change, and then they further developed a non-dominated sorting-based GA for multi-objectives optimization of the material composition of a three-parameter FG beam, in which the TI and TD material properties are considered. All of the results in the tables and figures are obtained using the self-developed MATLAB programs without any commercial toolbox.
In the thermal buckling analysis, the results show that the TD material properties must be considered in TD physical problems, in which for a moderately thick FG beam, the deviations in the critical temperature change parameters obtained using the TD material properties and the TI material properties are as high as approximately 36.5%, and the analysis based on the TI material materials is an unsafe analysis, due to the fact that it always overestimates the actual critical temperature change parameters.
In the optimal design cases, the results show the self-developed non-dominated sorting-based GA converges rapidly, where the Poreto-optimal solutions can be yielded after the 20th generation. The non-dominated sorting-based GA can also be extended to other optimal designs of FG beams with multiple objective functions.  Informed Consent Statement: This study did not involve humans and animals. Data Availability Statement: All of the data for Figures 1-9 can be obtained at https://drive.google. com/drive/folders/1J2KY29JdaCAhlpuGJdLSz0fb3kU_Z-HP?usp=sharing (accessed on 30 March 2021).

Conflicts of Interest:
The authors declare no conflict of interest.