A Multi-Level Optimization Method for Elastic Constants Identiﬁcation of Composite Laminates

: A new and e ﬀ ective elastic constants identiﬁcation technique is presented to extract the elastic constants of a composite laminate subjected to uniaxial tensile testing. The proposed technique consists of a new multi-level optimization method that can solve di ﬀ erent types of minimization problems, including the extraction of material constants of composite laminates from given strains. In the identiﬁcation process, the optimization problem is solved by using a stochastic multi-start dynamic search minimization algorithm at the ﬁrst level in order to obtain the statistics of the quasi-optimal design variables for a set of randomly generated starting points. The statistics of the quasi-optimal elastic constants obtained at this level are used to determine the reduced feasible region in order to formulate the second-level optimization problem. The second-level optimization problem is then solved using the particle swarm algorithm in order to obtain the statistics of the new quasi-optimal elastic constants. The iteration process between the ﬁrst and second levels of optimization continues until the standard deviations of the quasi-optimal design variables at any level of optimization are less than the prescribed values. The proposed multi-level optimization method, as well as several existing global optimization algorithms, is used to solve a number of well-known mathematical minimization problems to verify the accuracy of the method. For the adopted numerical examples, it has been shown that the proposed method is more e ﬃ cient and e ﬀ ective than the adopted global minimization algorithms to produce the exact solutions. The proposed method is then applied to identify four elastic constants of a [0 ◦ / ± 45 ◦ ] s composite laminate using three strains in 0 ◦ , 45 ◦ , and 90 ◦ directions, respectively, of the composite laminate subjected to uniaxial testing. For comparison purposes, several existing global minimization techniques are also used to solve the elastic constants identiﬁcation problem. Again, it has been shown that the proposed method is capable of producing more accurate results than the adopted available methods. Finally, experimental data are used to demonstrate the applications of the proposed method.


Introduction
As is well known, the determination of the true elastic constants is an essential step for design, fabrication, and the structural health monitoring of structures. Therefore, the development of efficient and effective methods for determining the actual elastic constants of materials and structures has become an important topic of research. However, for engineering applications, having simple testing procedures that can produce good estimations of the elastic constants of composite materials in an efficient way is always desirable. In recent years, measured mechanical characteristics, together with optimization techniques such as the genetic algorithm, particle swarm optimization method, stochastic laminates can produce good results, in terms of time, it is still expensive to fabricate the laminates and perform the tests. Therefore, as far as engineering applications are concerned, it is worthwhile to investigate whether it is possible to use only one composite laminate in the simple uniaxial tensile testing to identify the elastic constants of the laminate. It is noted that if a non-angle-ply composite laminate is used in the identification process, the search for accurate elastic constants may become very difficult and the previously proposed optimization techniques for elastic constants identification using one or two angle-ply laminates may be unable to produce good results in an efficient and effective way. Therefore, a new and effective optimization technique should be developed to extract elastic constants of a non-angle-ply composite laminate.
In this paper, a simple yet effective multi-level optimization method for elastic constant identification of composite materials is presented. The solutions of a number of well-known minimization problems using the proposed multi-level optimization method, as well as several other existing global minimization algorithms, are first given to demonstrate the capability of the proposed method to search for the global minima. The elastic constants identification of a [0 • /±45 • ] s composite laminate is then performed to illustrate the applications of the proposed method. The accuracy and feasibility of the proposed method to identify elastic constants of different composite materials are studied by way of several numerical and experimental examples.

Multi-Level Optimization Method
The following constrained minimization problem should be considered, which is stated in terms of N design variables x = = [x 1 , x 2 , . . . ., x N ], the objective function F(x), and a set of initially prescribed side constraints g(x) = [g 1 (x), g 2 (x), . . . ., g 2N (x)] T ≤ 0, with the superscript T denoting the transpose of the row vector.
Minimize F(x) Subject to g k (x) ≤ 0; k = 1, 2, . . . , 2N ( In general, a number of local minima in the feasible region of the above minimization problem may exist, especially in the case of elastic constants identification of composite materials. Herein, a multi-level optimization method (MLOM) in which several level-wise optimization problems were solved sequentially for determining the global minimum was proposed to deal with the above constraint minimization problem. In the proposed method, for any level of optimization, a set of starting points, randomly generated from a uniform distribution in the level-wise feasible region, were used to determine the quasi-global minima when solving the level-wise minimization problem. The statistics (means and standard deviations) of the quasi-optimal design variables associated with the quasi-global minima were determined to see if the true global minimum was attained, and if not, whether it was necessary to construct the new feasible region of the next level of optimization. The search for the true global minimum at different levels was terminated when the standard deviations of the quasi-optimal design variables were less than the prescribed value. Without the loss of generality, the commonly used global minimization algorithms such as the multi-start global search technique (MGST) [20] and particle swarm optimization technique (PSOT) [21] were adopted, respectively, at two consecutive levels to determine the quasi-global minima. Hence, in the proposed MLOM, the MGST was used to solve the first level optimization problem. The means and standard deviations of the quasi-optimal design variables obtained at this level were used to produce a contracted/reduced level-wise feasible region to establish the second level optimization problem. Following this line, the MGST and PSOT took turns to solve the subsequent levels of optimization problems. The advantage of the proposed solution strategy was that the adoption of a series of contracted level-wise feasible regions at different levels enhanced the chance of attaining the true global minimum. Herein, at the first-level of optimization, the constrained minimization problem of Equation (1) was converted to the following unconstrained minimization problem by introducing the general augmented Lagrangian [22].
where the superscript U and subscript L denote the upper and lower bounds, respectively; α j , η j , and r p are multipliers; and max [*,*] takes on the maximum value of the numbers in the bracket. The updated formulas for the multipliers α j , η j , and r p are where the superscript n denotes the iteration number; γ 0 is a constant; and r max p is the maximum value of r p . Following the guidelines given in the literature [22], the parameters µ 0 j , η 0 j , r 0 p , γ 0 , and r max p were chosen as The above unconstrained minimization problem was solved using the MGST. First, a set of starting points were randomly generated and for each starting point, the MGST was used to search for the lowest local minimum along the search path. The adopted search technique was formulated on the basis of the energy conservation of a travelling point mass, which tends to roll down from a hilltop (starting point) towards the valleys (local minima) along a chosen trajectory in the feasible region. The search did not stop until the rolling point mass reached a position/height with potential energy equal to that of the starting point. Therefore, a number of local minima could be attained along the search trajectory and the comparison among these local minima produced the lowest local minimum, which was termed a quasi-global minimum. In this way, a set of quasi-global minima could be obtained for the adopted starting points. Through a series of numerical tests, it was found that the quasi-global minima were generally in the vicinity of the global minimum. The statistics, namely, the means, µ i , and standard deviations, σ i , of the design variables associated with the quasi-global minima were calculated to determine whether the global minimum was attained, and if not, the current feasible region was updated. Herein, the means of the design variables associated with the quasi-global minima were defined as the quasi-optimal design variables at the level of optimization under consideration, while the standard deviations were used to measure the dispersions of the quasi-optimal design variables. It was assumed that, based on the limit theorem, when the number of starting points was large, the distributions of the quasi-optimal design variables were normal, and the expected values of the quasi-optimal design variables were the true design variables. It is noted that the quasi-global minimum was treated as the global minimum when the standard deviations of the quasi-true design variables are less than the small-prescribed value σ ζ . Otherwise, the quasi-optimal design variables and their standard deviations were used to construct the contracted feasible range, which became the level-wise feasible range of the second-level optimization problem. The contracted level-wise feasible ranges of the design variables for the next level were obtained as where q i are constants. Based on a series of numerical analyses, as will be described in the following, the choice of [1 ≤ q i ≤ 2] was capable of producing appropriately reduced feasible regions. Once the (reduced) level-wise feasible ranges were available, the second-level optimization problem could be stated as where g k * (x) ≤ 0 are the contracted side constraints derived from Equation (6).
A detailed description of the solution of a constrained optimization problem using PSOT can be found in the literature [23,24]. Herein, a brief introduction of the PSOT used to search for the global minimum is given. The above constrained optimization problem is first converted to the unconstrained minimization problem, as stated in Equation (2) and the M sets of K particles (starting points) that were randomly generated in the feasible region of this level were used in the PSOT to search for the global minimum. There were M runs to be performed. For each run, a set of K particles was used in the process of searching for the global minimum. The search involved the determination of the new positions of the particles via an iterative approach. For illustration, the search process of a particle can be used to show the idea of updating the position of the particle, as shown schematically in Figure 1. The movement of the particle in the iterative process of global minimum search is governed by the following equations: where V i and X i illustrate the speed and position of particle i, respectively; t is the time step; ω is inertia weight; c 1 is the cognitive coefficient; c 2 is the social coefficient; r 1 and r 2 are the random coefficients for cognitive and social coefficients, respectively; pBest is the individual optimal solution of the particle at time step t; gBest is the current global optimal solution of all the particles in the chosen set before the time step t + 1. In each iteration step, the position of the particle at time step t + 1 was derived from the information of the position and speed of the particle at time step t and two reference particle positions, i.e., one was the optimal solution of the particle itself (individual optimal solution, pBest) at time step t and the other was the temporary global optimal solution of all the particles (group optimal solution, gBest) before time step t + 1. It is noted that the use of a large inertia weight induces large particle motion inertia, of which the consequence is that the particle swarm will have a strong global exploration, i.e., a wide range of global search. On the contrary, the use of a small inertia weight will restrict the particle motion inertia, which makes the particle swarm have a strong local search exploration. Since the size of the reduced level-wise feasible region adopted for this level of optimization has become small, a small inertia weight was used to search for the quasi-global minimum of the K particles. It is noted that during the global minimum search, if the pBest of a particular particle, which is not the current group global minimum (gBest), repeatedly takes on the same local minimum for N p times, a new particle will be generated to replace that particle. The quasi-global minima of the M sets of particles will lead to M sets of optimal design variables, which can be used to determine statistics such as the means and standard deviations of the optimal design variables. Again, the means of the M sets of optimal design variables are termed as quasi-optimal design variables and a new reduced feasible region was attained by following the rule in Equation (6). It is noted that the goal was achieved and the means of the quasi-optimal design variables were treated as the true design variables when the standard deviations of the quasi-true design variables were less than the prescribed value σ ζ . Otherwise, the search for the global minimum continued by returning back to the first level of optimization with the use of the new contracted feasible region. Therefore, the iteration process between the first and second levels of optimization continued until all the standard deviations of the quasi-optimal design variables at any level of optimization were less than the prescribed value σ ζ .
www.mdpi.com/journal/applsci standard deviations of the quasi-true design variables were less than the prescribed value σζ. Otherwise, the search for the global minimum continued by returning back to the first level of optimization with the use of the new contracted feasible region. Therefore, the iteration process between the first and second levels of optimization continued until all the standard deviations of the quasi-optimal design variables at any level of optimization were less than the prescribed value σζ. Figure 1. The schematic diagram of the particle swarm optimization technique.

Numerical Examples of Global Minimum Determination
The proposed MLOM was first used to find the global minima of the following well-known functions. (9).
The function f in Equation (9) or (10) that was minimized can be expressed in the form of Equation (1). The proposed multi-level optimization technique, as described in Section 2, was used to search for the global minimum of the function. At the first level where MGST was applied to solve the unconstrained minimization problem, as stated in Equation (2), ten starting points were randomly generated from a uniform distribution to determine the quasi-global minima. For example, the contracted feasible region, adopted for the next level, was determined as [μi ± qiσi], The schematic diagram of the particle swarm optimization technique.

Numerical Examples of Global Minimum Determination
The proposed MLOM was first used to find the global minima of the following well-known functions.
for − 15 x −5, −3 y 3 or for − 100 x 100, −100 y 100 (10) The function f in Equation (9) or (10) that was minimized can be expressed in the form of Equation (1). The proposed multi-level optimization technique, as described in Section 2, was used to search for the global minimum of the function. At the first level where MGST was applied to solve the unconstrained minimization problem, as stated in Equation (2), ten starting points were randomly generated from a uniform distribution to determine the quasi-global minima. For example, the contracted feasible region, adopted for the next level, was determined as [µ i ± q i σ i ], where q i = 1. At the second level, where PSOT was applied to solve the optimization problem, five sets, of which each set contained 10 randomly generated starting points, were used to search for the quasi-global minima. The contracted feasible region, adopted for the next level, was also determined as [µ i ± σ i ]. The iteration between the two levels continued until the standard deviations of the quasi-true design variables were less than the prescribed value (σ ζ = 10 −7 ). Herein, three levels were used to search for the global minima of the functions. Using the Bukin function (I) as an example, the initial and reduced feasible regions at different levels of optimization are shown in Figure 2. It is noted that, at any level, the reduced feasible region was significantly smaller than the initial one. www.mdpi.com/journal/applsci sets, of which each set contained 10 randomly generated starting points, were used to search for the quasi-global minima. The contracted feasible region, adopted for the next level, was also determined as [μi ± σi]. The iteration between the two levels continued until the standard deviations of the quasi-true design variables were less than the prescribed value (σζ = 10 −7 ). Herein, three levels were used to search for the global minima of the functions. Using the Bukin function (I) as an example, the initial and reduced feasible regions at different levels of optimization are shown in Figure 2. It is noted that, at any level, the reduced feasible region was significantly smaller than the initial one.  The solutions of the functions, including the statistics of the identified design variables, are listed in Table 1. It is noted that the present method can very accurately predict results with less than or equal to 0.01% errors. For comparison purposes, other global minimization techniques, for instance, genetic algorithm (GA) and particle swarm (PS), available in MATLAB and the Statistical Multi-start Global Optimization Method (MSGO) [25], were also used directly to search for the The solutions of the functions, including the statistics of the identified design variables, are listed in Table 1. It is noted that the present method can very accurately predict results with less than or equal to 0.01% errors. For comparison purposes, other global minimization techniques, for instance, genetic algorithm (GA) and particle swarm (PS), available in MATLAB and the Statistical Multi-start Global Optimization Method (MSGO) [25], were also used directly to search for the global minima of the functions. The GA and PS are available in the global optimization toolbox of MATLAB [26]. Regarding GA, ten sets of randomly generated starting points were used to find the global minimum. For each set, the number of starting points (population size) was 10,000, the maximum number of generations was 10,000, the crossover fraction was 0.8, and the mutation rate was 10 −5 . Here, the true global minimum was treated as the lowest value of the quasi-global minima, which was identified from the ten sets of randomly generated starting points. Regarding PS, ten sets, of which each set had 100 randomly generated particles, were chosen to find the global minimum with the use of the following coefficients: The coefficient of inertia 0.01, the cognitive coefficient 0.1, and the social coefficient 0.01. Again, the true global minimum was treated as the lowest value of the quasi-global minima identified from the ten sets. Regarding MSGO, 100 starting points were randomly generated to search for the global minimum with the acceptance probability equal to 0.98. The identified design variables predicted using different methods are listed in Table 2 for comparison. It is noted that the overall performance of the present method was superior to the other methods. For instance, among the adopted optimization methods, only PS and MLOM could produce the exact values for the design variables of the Bukin function (I) with small bounds of the design variables. However, for larger bounds, e.g., the Bukin function (II), only the present MLOM could produce the exact solution, while the prediction using PS had errors as large as 5%. Furthermore, it is noted that the total number of starting points generated in the proposed method was much less than those of the selected global minimization methods. This revealed the fact that the proposed method was efficient and capable of producing accurate results. Table 1. Global minimum of function predicted using the proposed multi-level optimization method (MLOM).

Function
Level

Elastic Constants Identification of Composite Laminate
To demonstrate one of its applications, the proposed MLOM was used to identify the elastic constants of symmetric composite laminates, in which the composite layers had the same material properties, but may have had different fiber orientations. Herein, strains on the surface of the composite laminate were used to determine the four elastic constants (E 1 , E 2 , G 12 , and ν 12 ) of the laminate.

Strain Analysis of Symmetric Composite Laminate under Axial Load
The symmetric laminate of size a(length) × b(width) × h(thickness) subjected to in-plane axial stress resultant N x in the x-direction should be considered, as shown in Figure 3. The relations between the axial stress resultant and strains expressed in matrix form are [27] where ε x , ε y , and γ xy are axial, lateral, and engineering shear strains, respectively; A xx , A yy , A xy , and A ss are in-plane axial, lateral, axial-lateral coupling, and shear stiffness coefficients, respectively. The total force in the x-direction is defined as F x = bN x .

Elastic Constants Identification of Composite Laminate
To demonstrate one of its applications, the proposed MLOM was used to identify the elastic constants of symmetric composite laminates, in which the composite layers had the same material properties, but may have had different fiber orientations. Herein, strains on the surface of the composite laminate were used to determine the four elastic constants (E1, E2, G12, and ν12) of the laminate.

Strain Analysis of Symmetric Composite Laminate under Axial Load
The symmetric laminate of size a(length) × b(width) × h(thickness) subjected to in-plane axial stress resultant Nx in the x-direction should be considered, as shown in Figure 3. The relations between the axial stress resultant and strains expressed in matrix form are [27] 0 00 0 0 0    The laminate in-plane stiffness coefficients are where h is laminate thickness and Q with where E 1 and E 2 are Young's moduli in fiber and transverse directions, respectively; ν ij is Poisson's ratio for transverse strain in the j-direction when stressed in the i-direction; and G 12 is shear modulus in the 1-2 plane. The relations between Q ij and Q ij can be found in the literature. (11) gives

The inversion of Equation
a xx a xy 0 a xy a yy 0 0 0 a ss where a ij (i, j = x, y, s) are the in-plane compliance coefficients, which are functions of the elastic constants (E 1 , E 2 , G 12 , and ν 12 ) and the layer fiber angles of the laminate. Hence, from the above equation, the axial and lateral strains can be obtained, respectively, as and ε y = a xy N x (17) It is noted that the strains (ε x and ε y ) can be determined directly from Equations (16) and (17), provided that the elastic constants of the composite laminate are available. However, the strains of a composite laminate can generally be determined experimentally without knowing the values of the elastic constants. It is a premise that through an optimization approach, the elastic constants of a laminate may be extracted from the strains of the laminate, even though the number of strains was less than the number of elastic constants. Herein, as an example, the extraction of the elastic constants (E 1 , E 2 , ν 12 , and G 12 ) from given strains of a symmetric [0 • /±45 • ] s laminate will be achieved via an optimization approach. Traditionally, the given or measured axial and lateral strains (ε * x and ε * y ), together with the theoretically predicted (ε x and ε y ), can be used to construct the objective function to formulate the optimization based elastic constants identification problem. For instance, the commonly used least square type objective function was established as the weighted sum of (ε * x − ε x ) 2 and (ε * y − ε y ) 2 . Through a detailed numerical study, it has been shown that, if the given or measured strains are exact, i.e., contain no noise, the use of the above least square type objective function, comprised of two given strains in the optimization based elastic constants identification, can produce excellent results. However, for practical consideration, uncertainties involved in the measurement of strains during testing are inevitable. If there is some noise in the measured strains (ε * x and ε * y ), the identified elastic constants might deviate from their true values. Therefore, it was desirable to adopt a strain constraint to lower the dominance of the measured axial or lateral strain and regulate the variations of the theoretical axial and lateral strains so that the noise effects could be minimized. Before proceeding to the elastic constants identification, it is noted that, for a 0 • or 90 • layer with given elastic constants, the theoretically predicted normal strain in 45 • direction, ε 45 , could be calculated directly from the theoretical axial and lateral strains of the layer as Herein, the given or measured normal strain in 45 • direction, ε * 45 , and the calculated ε 45 , was used to establish the strain constraint. The role of the strain constraint was to regulate the motions of the design variables so that the theoretical normal strain in 45 • direction, calculated from Equation (18), was forced to have a good match with the given normal strain in 45 • direction. Therefore, the given strains (ε * x , ε * y and ε * 45 ), together with the theoretical strains (ε x , ε y and ε 45 ), were used in the proposed MLOM to identify the elastic constants.

Formulation of Elastic Constants Identification
The minimization problem of elastic constants identification is expressed as where e(x), containing three strain difference terms, is the objective function measuring the differences between the predicted and given (measured) strains; x = [E 1 , E 2 , G 12 , ν 12 ] the elastic constants with , and x 4 = ν 12 ; x L i and x U i are, respectively, the lower and upper bounds of the elastic constants x i ; ξ i (i =1-3) are weighing factors used to make the strain components have appropriate contributions to the error function and avoid the occurrence of numerical underflow of the error function. In view of Equation (18), when ε 45 was expressed in terms of ε x and ε y , the objective function in Equation (19), comprising the coupling terms of ε x and ε y , was no longer a least square type objective function. The proposed MLOM, as described in Section 2, was then used to solve the above minimization problem.

Numerical Examples
As an example, the strains of a Graphite/epoxy [0 • /±45 • ] s composite laminate were used in the proposed multi-level optimization method to identify the elastic constants of the constituent composite material. The actual elastic constants of the composite material were set as E 1 = 138.5 GPa, E 2 = 8.1 GPa, G 12 = 5.2 GPa, and ν 12 = 0.31. When subjected to axial force F x = 0.5 kN, the actual strains of the laminate computed from Equations (16)- (18) were ε x (µε) = 210.49, ε y (µε) = −151.81, and ε 45 (µε) = 29.34. Herein, the actual strains were treated as the "given" strains, which were used to identify the true elastic constants. To cover different kinds of fiber-reinforced composite materials, the upper and lower bounds of the elastic constants were initially chosen as At the first level of optimization, the values of the weighting factors ξ i (i = 1-3) in Equation (19) could be chosen arbitrarily so long as these values were positive and the weighting factors of the strain difference terms in Equation (19), containing the axial strain, had larger values than that of the strain difference term that only contained the lateral strain. Herein, based on engineering judgment, the values of the weighting factors ξ i (i = 1-3) were set to be 10 6 , 10 4 , and 10 6 , respectively. This choice could have led to more accurate prediction of E 1 . On the contrary, at the second level of optimization, the values of the weighting factors ξ i (i = 1-3) in Equation (19) were set to be 10 4 , 10 6 , and 10 4 , respectively. At this level, the larger weight for the lateral strain may have led to more accurate predictions of E 2 , while only a small adjustment of E 1 might occur. The procedure of the elastic constants identification of the graphite/epoxy laminate was first described in detail to illustrate the application of the MLOM. At the first level, the constrained minimization problem of Equation (19) was converted to an unconstrained minimization one, as stated in Equation (2). The unconstrained minimization problem was then solved using the MGST, in which ten starting points were randomly generated in the initial feasible range to determine the quasi-global minima and their associated quasi-optimal design variables. The starting points and the statistics of the quasi-optimal design variables are listed in Table 3. It is noted that the means of the quasi-optimal design variables were very close to the actual values, with errors less than or equal to 3.19%. The contracted feasible region adopted for the next level was determined as [µ i ± σ i ]. At the second level where PSOT was applied to solve the optimization problem, five sets, of which each set contained ten randomly generated starting points, were used to search for the quasi-global minima. The statistics of the quasi-optimal design variables are listed in Table 4. The comparison between the previous and current levels shows that the means of the quasi-optimal elastic constants of the current level were much closer to the actual values, with errors less than or equal to 1.51%. The contracted feasible region adopted for the next level was also determined as [µ i ± σ i ]. The iteration between the two levels continued until the standard deviations of the quasi-true design variables were less than a prescribed value (σ ζ = 10 −3 ). Finally, the solution of the optimization problem converges at the third level. The identified elastic constants associated with the global minimum are listed in Table 5 in comparison to the actual elastic constants. It is noted that the differences between the actual and identified elastic constants were less than or equal to 0.97%. To demonstrate the advantage of the proposed method, the commonly used global minimization methods were also used to identify the elastic constants of the graphite/epoxy [0 • /±45 • ] s composite laminate. The values of the parameters used in the selected global minimization methods have been given in Section 2. The elastic constants identified using the selected global minimization methods are listed in Table 6 in comparison to those identified using the proposed method. It is clear that the accuracy of the proposed method was much higher than the other methods. Again, it is noted that the total number of starting points generated in the proposed method was much less than those of the adopted popular global minimization methods. This revealed the fact that the proposed method was efficient and capable of producing accurate results. To further demonstrate the capability of the proposed MLOM, the elastic constants of different composite materials were also extracted from their given strains. The actual elastic constants, together with the given strains used for elastic constants identification of two composite materials, namely, carbon/epoxy and glass/epoxy, respectively, are listed in Table 7. The solutions of the first and second levels of optimization problems are listed in Table 8. Again, it is noted that at the second level, the mean of the quasi-optimal elastic constants were very close to the actual elastic constants with errors less than or equal to 0.09% for the carbon/epoxy and 1.5% for the glass/epoxy. At the third level, the identified elastic constants for the composite materials are listed in Table 9 in comparison to the actual ones. It is noted that the errors of the identified elastic constants were less than or equal to 0.06% and 0.99% for the carbon/epoxy and glass/epoxy, respectively. The small errors of the identified elastic constants revealed the capability of the proposed MLOM for elastic constants identification of different composite materials.

Experimental Study
The proposed method was then applied to the elastic constants of a symmetric graphite/epoxy [0 • /±45 • ] s laminate using three measured strains of the laminate. First, the material constants of the graphite/epoxy lamina were determined using three types of the standard specimens ([0 • ] 4 , [90 • ] 8 , and [±45 • ] s ), in accordance with the ASTM standards of D3039 and D3518 [28]. The average values and coefficients of variation (indicated in the parentheses) of the material constants are given as follows: In view of the small coefficients of variation (less than 2%), the average values in the above equation were then treated as the "actual" elastic constants of the composite material.
The dimensions of the symmetric graphite/epoxy [0 • /±45 • ] s laminate are shown in Figure 4. The average thickness of each lamina in the laminate is 0.269 mm. The graphite/epoxy [0 • /±45 • ] s laminates were subjected to uniaxial tensile testing, in which a tri-axial strain gauge produced by TML, Japan (TML FRA-6-11-1L, gage length 6mm and gauge factor 2.11 ± 1%) were used to measure the surface axial, lateral, and 45 • strains (ε * x , ε * y , and ε * 45 ) at the mid-span of the laminate. The uniaxial tensile test of the laminate was performed using the MTS 810, with test speeds less than 0.01 mm/s. The load-strain relation of the laminate was constructed to determine the strains for any given load and the elastic constants of the constituent material. The measured axial, lateral, and 45 • strains of the laminates for axial load F x = 0. 5  In view of the small coefficients of variation (less than 2%), the average values in the above equation were then treated as the "actual" elastic constants of the composite material.
The dimensions of the symmetric graphite/epoxy [0°/±45°]s laminate are shown in Figure 4. The average thickness of each lamina in the laminate is 0.269 mm. The graphite/epoxy [0°/±45°]s laminates were subjected to uniaxial tensile testing, in which a tri-axial strain gauge produced by TML, Japan (TML FRA-6-11-1L, gage length 6mm and gauge factor 2.11 ± 1%) were used to measure the surface axial, lateral, and 45° strains ( * , * , and 45 * ) at the mid-span of the laminate. The uniaxial tensile test of the laminate was performed using the MTS 810, with test speeds less than 0.01 mm/s. The load-strain relation of the laminate was constructed to determine the strains for any given load and the elastic constants of the constituent material. The measured axial, lateral, and 45° strains of the laminates for axial load Fx = 0.5 kN are ε * (με) = 208.57 (0.91%) ε * (με) = −149.89 (1.26%) ε 45 * (με) = 28.86 (1.64%) (22). For comparison, the "actual" elastic constants were used in Equations (16)-(18) to calculate the "actual" strains. The absolute percentage differences between the "actual" and measured strains, as shown in the parentheses, are given in the above equation. It is noted that the deviations of the measured strains from the "actual" ones are likely due to the existence of uncertainties in material properties, fabrication of specimens, strain measurement, etc. These uncertainties and their effects should be studied in detail in the future. The measured strains were then used in the proposed multi-level optimization method to identify the elastic constants. The identified elastic constants and the percentage differences between the "actual" and identified elastic constants are given as  For comparison, the "actual" elastic constants were used in Equations (16)-(18) to calculate the "actual" strains. The absolute percentage differences between the "actual" and measured strains, as shown in the parentheses, are given in the above equation. It is noted that the deviations of the measured strains from the "actual" ones are likely due to the existence of uncertainties in material properties, fabrication of specimens, strain measurement, etc. These uncertainties and their effects should be studied in detail in the future. The measured strains were then used in the proposed multi-level optimization method to identify the elastic constants. The identified elastic constants and the percentage differences between the "actual" and identified elastic constants are given as E 1 (GPa) = 140.10(1.15%) E 2 (GPa) = 8.093(0.086%) v 12 = 0.302(0.026%) G 12 (GPa) = 5.185(0.29%) (23) It is noted that, even though the accuracy of the measured strains was affected by uncertainties, the proposed multi-level method could still produce excellent results with absolute percentage differences less than or equal to 1.15%. Furthermore, it is worthy to note that, if the constraint term of (ε * 45 − ε 45 ) 2 is not included in the objective function of Equation (19), the percentage differences between the "actual" and identified elastic constants might exceed 20%. Therefore, the strain constraint term contributed by making the elastic constants refrain from having large deviations from the actual values. The excellent results obtained for the [0 • /±45 • ] s laminate might lead to the possibility of having a new testing procedure for elastic constants identification in the future. However, the uncertainties involved in laminate fabrication and strain measurement should be characterized and their effects on the elastic constants analyzed before this new testing procedure can be standardized.

Conclusions
An efficient and effective multi-level optimization method, consisting of several levels of optimization, was presented to solve different types of minimization problems. The proposed optimization method included two major parts. Regarding the first part, the initial minimization problem was tackled via the solutions of a series of level-wise minimization problems, in which the level-wise feasible region of any level was derived from that of the previous level via a statistical approach. As for the second part, the multi-start global search technique (MGST) and the particle swarm optimization technique (PSOT) are used, respectively, to solve two consecutive level-wise minimization problems. The locations of the global minima of several functions were identified using the proposed method, as well as several commonly used global minimization algorithms. The total number of starting points generated in the proposed method was much less than those generated in each of the selected minimization methods. It has been shown that the proposed method could produce more accurate results than the selected methods. To demonstrate its applications, the proposed optimization method has been applied to the elastic constants identification of composite laminates subjected to a uniaxial load. Three measured strains in 0 • , and 45 • , and 90 • directions have been used to establish the objective function. It has been shown that the measured normal strain in the 45 • direction can work as a constraint to regulate the variations of the theoretical axial and lateral strains, and in turn make the design variables move efficiently toward their optimal values. Therefore, more accurate results can be obtained if the measured 45 • strain is taken into account in establishing the objective function. The elastic constants of [0 • /±45 • ] s laminates made of different composite materials have been identified using the 0 • , 45 • , and 90 • normal strains on the top layer of each laminate. For comparison, the selected global minimization techniques have also been used to perform the elastic constants identification. Again, the total number of starting points generated in the proposed method was much less than those generated in each of the selected global minimization methods. It has been shown that the proposed method with the use of three levels and a few starting points randomly generated at each level can produce more accurate results than the selected global minimization techniques. The absolute errors of the elastic constants identified are less than or equal to 0.97%, 0.06%, and 0.99% for the graphite/epoxy, carbon/epoxy, and glass/epoxy, respectively. Finally, the tensile test of a graphite/epoxy [0 • /±45 • ] s laminate was performed. The proposed method was applied to identify the elastic constants of the laminate using three measured strains. In the face of existing uncertainties, excellent results with a percentage of differences equal to or less than 1.15% have been obtained to reveal the potential of the proposed method for practical applications.
Author Contributions: C.Y.H. did the mathematical programming, analyses, experimental work, and data collection. T.Y.K. set the goal, constructed the method, and wrote the paper.