Inverse identification of elastic properties of composite materials using hybrid GA-ACO-PSO algorithm

Abstract

The main emphasis of this paper is placed on the effectiveness of the proposed optimization method in material identification. The primary motivation of integrating GA, ACO and PSO is to minimize each other’s weaknesses and to promote respective strengths. In the proposed algorithm, the effect of random initialization of GA is subdued by passing the products of GA through the ACO and PSO operators to well organize the exploitative and exploratory search coverage. In return, GA improves the convergence rate and alleviates the strong dependency on the pheromone array in ACO as well as resolves the conflict arisen in identifying the trade-off parameter and further refine the exploitative search of PSO with the introduction of two-point standard mutation and one-point refined mutation. The proposed algorithm has been verified and applied in composite material identification with absolute percentage errors between measured and evaluated natural frequencies not more than 2%.

Keywords:

• Ant colony optimization
• composites
• genetic algorithm
• hybrid algorithm
• particle swarm optimization

AMS Subject classifications:

• 34K29
• 65K10

1. Introduction

Great interests have been focalized on developing non-destructive approaches due to the ease of implementation, rapid acquisition and cost-saving factor. This approach is usually referred to as experimental–numerical method as it involves the acquisition of experimental data, which are then used to numerically evaluate the material properties of a specimen. In vibration, experimental modal analysis (EMA) is conducted to study the dynamic behaviours of a structure, i.e. natural frequencies, mode shapes and damping properties. Based on either/all of these acquired experimental modal parameters, numerical algorithm is subsequently used to evaluate the elastic properties of the material. Destruction of specimen is not required in this approach. Due to its prominent advantages, this experimental–numerical identification method has been widely used in recent decades. In the context of numerical evaluation, it comprises forward methods and inverse methods. Elastic properties are used as inputs in forward method to evaluate modal parameters, whereby, in inverse methods, error function(s) is defined and optimized (minimization or maximization) involving the correlation difference between experimental and evaluated modal parameters. The accuracy of a non-destructive approach depends greatly on the accuracy of both forward methods as well as the inverse or optimization methods. Rayleigh’s method [1,2], Rayleigh–Ritz method [3,4] and finite element method [5–7], are amongst the most prevalently used forward methods due to numerous advantages. Nevertheless, the main drawback of Rayleigh’s method lies in its inability of providing higher modes of modal solutions, while for Rayleigh–Ritz method, the difficulty in predicting the appropriate shape functions for specific problems impedes its usage preference in tackling complex problems. Finite element method appears to be the solution to those mentioned limitations; however, if compared to Fourier method [8,9], its accuracy comes second. On the other hand, various inverse methods have been developed and utilized in material identification in the past as reviewed by Tam et al. [10]. Chakraborty and teammates [11] adopted the use of a basic GA to identify the geometric and material properties of isotropic and composite plates on account of its ability to perform simultaneous search for multiple solutions. Instead of opting basic GA, Han and Liu [12] later proposed the use of uniform μGA to determine the material properties of functionally graded material (FGM) plates. The difference between μGA and simple GA consists in the absence of mutation operator in μGA. It employs the technique of micropopulation which is able to evade local minima and exhibits higher rate of convergence, if compared to the conventional GA. To improve the quality of the solutions, Liu and co-workers [13] suggested the integrated use of μGA with a non-linear least squares method (LSM) to characterize the elastic properties of composite plates. The μGA operator was first used to search for better solutions close to the optima and these identified solutions were then used to initialize the non-linear LSM to locate the global optimum. Computational cost had been an issue, therefore, a modified μGA was proposed by Yang and Liu [14] to evaluate the material properties of a printed circuit board. An intergeneration projection searching technique (IP-μGA) was introduced to get rid of excessive number of functional evaluations, and thus, improving the convergence rate. Furthermore, hybrid real parameter GA was proposed by Hwang and colleagues [15,16], employing the concept of combining a simulated annealing mutation process and adaptive mechanisms with a real parameter GA to evaluate the elastic properties of composite plates. It showed excellent accuracy and practicality in composite material applications. Apart from GA, Majumdar and co-workers [17] utilized ACO to study damage of truss structures based on the changes in natural frequencies. Conventional ACO was opted in this study for its simplicity and effectiveness in solving such problem. Jalil et al. [18] developed a modified version of ant colony optimization (ACO) to construct a model representing a flexible plate structure. ACO with roulette wheel selection was proposed and the results revealed that the proposed ACO performed better, in terms of accuracy and convergence rate compared to previous version of ACO. Besides ACO, Samir [19] and Sankar [20] hybridized response surface methodology (RSM) approach with particle swarm optimization (PSO) method to improve the accuracy in identifying elastic properties of composite plates. Hybrid RSM-PSO imparted the most accurate results with the smallest errors amongst the different methods, while the convergence rate was reasonably good but not as good as that of solely RSM.

With regard to the reviews of hybrid algorithm that involves GA, ACO and/or PSO, Sheikhalishahi and co-authors [21] developed a hybrid GA-PSO approach to tackle reliability redundancy allocation problem. In the proposed approach, GA was used for more global search, while PSO was utilized for more precise search. However, only selected individuals were passed through PSO operator to reduce computational time. Yu and teammates [22] adopted the use of hybrid GA-PSO algorithm for process planning scheduling in manufacturing system. In the suggested algorithm, GA was assigned for the traversal of beginning process plans, followed by the use of PSO in choosing suitable machine for each process as well as figuring the initiation time of every process. Unlike the previous literature, PSO was applied on the entire populations. Lim and Ponnambalam [23] employed an integrated GA-PSO algorithm for the concurrent design of cellular manufacturing system. The novelty of this research can be seen in the simultaneous application of GA and PSO in making cell formation (CF) and grouping layout (GL) decisions, in which, GA was assigned for CF and PSO was designated for GL. Apart from GA-PSO, Fidanova and colleagues [24] suggested the use of hybrid GA-ACO algorithm to determine the parameters of a fermentation process. GA was initially used to search potential solutions in a coarse manner. Then, ACO was assigned to search for better solution in a more precise manner utilizing the information collected by GA. Furthermore, a GA-ACO-Local Search Hybrid algorithm was developed by Xu and collaborators [25] to solve quadratic assignment problem. The probability of GA and ACO taking place was set equal (0.5) and local search was then used to further process the product of either GA or ACO to locally search for better solutions. Besides GA-ACO, Teo and Ponnambalam [26] designed a combined ACO-PSO algorithm to tackle single row layout problem. In this hybrid algorithm, ACO was utilized initially to roughly search the feasible solutions, while PSO was then used to exploit the information garnered by ACO to search for better solutions at a more thorough pace. Gigras and Gupta [27] adopted the use of a hybrid ACO-PSO algorithm for robotic path planning. The concept of pheromone intensity in ACO and the concept of velocity in PSO were utilized collaboratively to identify the shortest path to improve the convergence rate.

Based on the aforementioned literatures, genetic algorithm (GA) is found to be one of the most renowned algorithms used in material identification and in industrial applications due to its key advantages. The salient advantage consists in its ability to return multiple solutions for further consideration. This is certainly important and useful when dealing with real-life problems containing multiple solutions. Furthermore, GA shows great potential in solving any problems that involve any types of objective functions. In the aspect of practicality, this will definitely be beneficial to those problems-prone industries, where, the presence of nonlinearity is inevitable. From a negative point of view, however, GA consumes much computational efforts when solving relatively complex problems and poses the possibility of falling into local minima. Pertaining to the strength of ACO, the concept of pheromone accumulation in ACO leads to more directive search such that the rate of opting repeated solutions are reduced over iterations, thus, this might reduce repeated evaluations as well as the overall number of evaluations. However, the main drawback consists in the strong dependency of the algorithm on the array of pheromone intensity that might lead to biased convergence. It might also consume more computational efforts and time when coping with complex problems. Furthermore, from the previous reviews, it can be observed that particle swarm optimization (PSO) is preferably adopted for collaborative application primarily because of its ease of implementation and high convergence rate. Nevertheless, its potential is only discernibly visible when the search region is limited to a certain range; else, the search might focalize on surrounding regions of a local solution instead of a global solution. Difficulty in determining the appropriate trade-off value that decides the distribution between exploratory search and exploitative search can be one of the challenges as well. Biased search will occur if the trade-off value is not properly selected, thus, might affect the outcomes. Hybrid algorithms involving the integration of GA-PSO, GA-ACO and ACO-PSO have been exhaustively developed and applied in various fields since decades ago for their great advantages. In most cases, GA is used for diversification and the other is utilized for intensification. Meanwhile, for the hybridization between ACO and PSO, ACO is normally used for exploration and PSO is assigned for exploitation. Although extensive efforts have been devoted to developing hybrids involving two basic algorithms, the application of hybrids encompassing more than two basic algorithms can hardly be found in material identification and thus, the collaboration between these three algorithms can be a potential approach.

Therefore, in this paper, the use of forward Fourier method is adopted for its proven prominent accuracy and a new hybrid algorithm, known as, GA-ACO-PSO algorithm is proposed and presented to identify the elastic properties of composite plates with different boundary conditions. The main emphasis of this paper is placed on the effectiveness of the inverse or optimization method in material identification. The primary motivation of integrating GA, ACO and PSO is to minimize respective weaknesses and in the meantime, to promote respective strengths and therefore, converting into a better hybrid algorithm. As mentioned, GA has the tendency of falling into local minima due to unorganized and random initialization. In the proposed algorithm, instead of eliminating the random initialization of GA, the effect of random initialization is subdued by passing the products of GA through the ACO and PSO operators, whereby, in the composition of the proposed algorithm, ACO is assigned for diversification, while PSO is designated for intensification. The products of GA will be further processed by ACO and PSO to well distribute and organize the search coverage and thus, neutralize the biased effect of random initialization. Logically, this will enhance the searching operation. Moreover, the presence of ACO operator with the concept of opting unvisited routes and the practice of pheromone accumulation can as well be the solution to counter the shortcoming of GA, where ACO is able to evade the occurrence of repeated and redundant tours. The concept of combined discrete and continuous probability distribution scheme of ACO algorithm is also employed to specifically assist and organize the exploratory search strategy in GA. As it is known that ACO has often been used for solving discrete problems, the suggested feature enables ACO to be used in solving continuous problems, whereby, the search space is initially divided discretely into several predefined search paths, followed by random and continuous search within each path. Furthermore, since GA exhibits inferiority in exploitation when comparing to that of PSO, PSO can thus be potentially utilized in only exploitative search to aid GA, where the search is brought into focus around potential solutions to identify the exact optimum in a precise manner. In return, GA will further enhance the convergence rate and reduce the strong dependency on the array of pheromone in ACO as well as help to overcome the conflict found in identifying the trade-off parameter and further refine the exploitative search strategy of PSO with the introduction of a different yet simple mechanism of crossover and mutation operators. Two types of mutation operators are suggested in this algorithm, namely, two-point standard mutation and one-point refined mutation. In early iterations, standard mutation will be utilized collaboratively with the concept of unrepeated tour in ACO to evade local entrapment, while refined mutation will be used in later iterations to supplement PSO intensification. The applicability of the proposed algorithm in material identification has been verified using examples taken from [28]. In the study, the plate dynamic model is constructed based on the use of Fourier series. Fourier method is used as the Forward method in this study for its excellent accuracy as proven by Khov [28]. The remaining paper is arranged as follows. Section 2 introduces the formulations of the plate dynamic model with general boundary conditions. Section 3 explains the general optimization objective function used in the study. Section 4 briefly introduces the fundamental concepts of conventional algorithms. Section 5 elaborates the proposed hybrid algorithm in details. Section 6 demonstrates the verification and application of the proposed algorithm in determining elastic properties and discussions are made. Section 7 concludes the paper based on the verified results and discussions.

2. Plate vibration model

In the present study, both shear deformation and rotary inertia are not taken into consideration as the span of the plate involved is much larger than its thickness (significantly large span-to-thickness ratio). As proven by Numayr et al. [29], a significant effect is only visible when shear deformation is considered in the computation of composite plate with small span-to-thickness ratio, while the inclusion of rotary inertia only shows a great influence when thick plates are involved. A plate is considered thick if the span-to-thickness ratio is less than 20 [30]. In the present study, the span-to-thickness ratios of the plates involved are more than 20 (about 100) and therefore, both shear deformation and rotary inertia are not considered in the present research to avoid unnecessary and complicated computation. The equation of motion for an orthotropic plate can be expressed as:(1) $$D_{11}\frac{{\mathit{\partial }}^{4}W}{\mathit{\partial }{x}^4}+4D_{16}\frac{{\mathit{\partial }}^{4}W}{\mathit{\partial }{x}^3\mathit{\partial }y}+2\left(D_{12}+2D_{66}\right)\frac{{\mathit{\partial }}^{4}W}{\mathit{\partial }{x}^2\mathit{\partial }{y}^2}+4D_{26}\frac{{\mathit{\partial }}^{4}W}{\mathit{\partial }x\mathit{\partial }{y}^2}+D_{22}\frac{{\mathit{\partial }}^{4}W}{\mathit{\partial }{y}^4}-\mathit{\rho }h{\mathit{\omega }}^{2}W(x,y)=q(\mathit{\omega },x,y)$$(1)

where ρ is the mass density of the plate, and ω is the frequency of the harmonic motion. (2) $$W\left(x,y\right)=\sum _{m=0}^{\infty }\sum _{n=0}^{\infty }{A}_{\mathit{mn}}cos\left({\mathit{\lambda }}_{\mathit{am}}x\right)cos\left({\mathit{\lambda }}_{\mathit{bn}}y\right)+\sum _{l=1}^4\left({}_b^l\left(y\right)\sum _{m=0}^{\infty }{c}_m^{l}cos\left({\mathit{\lambda }}_{\mathit{am}}x\right){+}_a^l\left(x\right)\sum _{n=0}^{\infty }{d}_n^{l}cos\left({\mathit{\lambda }}_{\mathit{bn}}y\right)\right)$$(2)

Equation (1(2) $$W\left(x,y\right)=\sum _{m=0}^{\infty }\sum _{n=0}^{\infty }{A}_{\mathit{mn}}cos\left({\mathit{\lambda }}_{\mathit{am}}x\right)cos\left({\mathit{\lambda }}_{\mathit{bn}}y\right)+\sum _{l=1}^4\left({}_b^l\left(y\right)\sum _{m=0}^{\infty }{c}_m^{l}cos\left({\mathit{\lambda }}_{\mathit{am}}x\right){+}_a^l\left(x\right)\sum _{n=0}^{\infty }{d}_n^{l}cos\left({\mathit{\lambda }}_{\mathit{bn}}y\right)\right)$$(2) ) describes the motion of an orthotropic plate, assuming that a harmonic excitation is initiated. A relationship between the stiffness rigidities, D _ij and the plate’s deflection displacement, $W\left(x,y\right)$ is established, in which the displacement deflection function is defined as a more robust form of Fourier series expansion, as shown in Equation (2). A system of equation can be derived by substituting Equation (2) into Equation (1(2) $$W\left(x,y\right)=\sum _{m=0}^{\infty }\sum _{n=0}^{\infty }{A}_{\mathit{mn}}cos\left({\mathit{\lambda }}_{\mathit{am}}x\right)cos\left({\mathit{\lambda }}_{\mathit{bn}}y\right)+\sum _{l=1}^4\left({}_b^l\left(y\right)\sum _{m=0}^{\infty }{c}_m^{l}cos\left({\mathit{\lambda }}_{\mathit{am}}x\right){+}_a^l\left(x\right)\sum _{n=0}^{\infty }{d}_n^{l}cos\left({\mathit{\lambda }}_{\mathit{bn}}y\right)\right)$$(2) ), and rearranged into Equation (3(3) $$\begin{array}{c}(D_{11}{\mathit{\lambda }}_{\mathit{am}}^4+D_{22}{\mathit{\lambda }}_{\mathit{bn}}^4+2\left(D_{12}+2D_{66}\right){\mathit{\lambda }}_{\mathit{am}}^2{\mathit{\lambda }}_{\mathit{bn}}^2)A_{\mathit{mn}}\hfill \\ +\sum _{l=1}^4(D_{11}{\mathit{\beta }}_n^l{\mathit{\lambda }}_{\mathit{am}}^4-2\left(D_{12}+2D_{66}\right){\mathit{\lambda }}_{\mathit{am}}^2{\overline{\mathit{\beta }}}_n^l+D_{22}{\overline{\mathit{\beta }}}_n^l)c_m^l\hfill \\ +\sum _{l=1}^4(D_{11}{\mathit{\alpha }}_m^l{\mathit{\lambda }}_{\mathit{bn}}^4-2\left(D_{12}+2D_{66}\right){\mathit{\lambda }}_{\mathit{bn}}^2{\overline{\mathit{\alpha }}}_m^l+D_{22}{\overline{\mathit{\alpha }}}_m^l)d_n^l\hfill \\ -\mathit{\rho }h{\mathit{\omega }}^2(A_{\mathit{mn}}+\sum _{l=1}^4({\mathit{\beta }}_n^l{c}_m^l+{\mathit{\alpha }}_m^l{d}_n^l))=q_{\mathit{mn}}\hfill \end{array}$$(3) ).(3) $$\begin{array}{c}(D_{11}{\mathit{\lambda }}_{\mathit{am}}^4+D_{22}{\mathit{\lambda }}_{\mathit{bn}}^4+2\left(D_{12}+2D_{66}\right){\mathit{\lambda }}_{\mathit{am}}^2{\mathit{\lambda }}_{\mathit{bn}}^2)A_{\mathit{mn}}\hfill \\ +\sum _{l=1}^4(D_{11}{\mathit{\beta }}_n^l{\mathit{\lambda }}_{\mathit{am}}^4-2\left(D_{12}+2D_{66}\right){\mathit{\lambda }}_{\mathit{am}}^2{\overline{\mathit{\beta }}}_n^l+D_{22}{\overline{\mathit{\beta }}}_n^l)c_m^l\hfill \\ +\sum _{l=1}^4(D_{11}{\mathit{\alpha }}_m^l{\mathit{\lambda }}_{\mathit{bn}}^4-2\left(D_{12}+2D_{66}\right){\mathit{\lambda }}_{\mathit{bn}}^2{\overline{\mathit{\alpha }}}_m^l+D_{22}{\overline{\mathit{\alpha }}}_m^l)d_n^l\hfill \\ -\mathit{\rho }h{\mathit{\omega }}^2(A_{\mathit{mn}}+\sum _{l=1}^4({\mathit{\beta }}_n^l{c}_m^l+{\mathit{\alpha }}_m^l{d}_n^l))=q_{\mathit{mn}}\hfill \end{array}$$(3)

In matrix form, it can be represented as shown in equation 4(4) $$\left(\tilde{K}A+BP)-\mathit{\rho }h{\mathit{\omega }}^2(\tilde{M}A+FP\right)=q$$(4) :(4) $$\left(\tilde{K}A+BP)-\mathit{\rho }h{\mathit{\omega }}^2(\tilde{M}A+FP\right)=q$$(4)

The natural frequencies and mode shapes of an orthotropic plate can then be obtained by solving a standard matrix characteristic Equation (5(5) $$\left(K-\mathit{\rho }h{\mathit{\omega }}^{2}M\right)A=q$$(5) ) with q = 0, in which the eigenvector A for a given eigenvalue contains the Fourier coefficients for the corresponding mode shapes. Using the boundary equations, HP = QA, Equation (4(4) $$\left(\tilde{K}A+BP)-\mathit{\rho }h{\mathit{\omega }}^2(\tilde{M}A+FP\right)=q$$(4) ) can be written in the form of(5) $$\left(K-\mathit{\rho }h{\mathit{\omega }}^{2}M\right)A=q$$(5)

where the stiffness matrix, $K=\tilde{K}+B{H}^{-1}Q$ and the mass matrix, $M=\tilde{M}+F{H}^{-1}Q$. The boundary conditions for this study are: F-F-F-F, C-F-F-F, C-C-F-C (F: free, C: clamped). The derivation and details of this forward method can be found in [28].

3. Material identification by the proposed hybrid GA-ACO-PSO algorithm

In the present study, the objective function can be defined in terms of discrepancy of resonant frequencies as follows:(6) $$\begin{array}{cc}\text{Minimise}\hfill & f\left(x\right)=\sum _{i=1}^N{\left(\frac{{\mathit{\omega }}_i-{\overline{\mathit{\omega }}}_i}{{\mathit{\omega }}_i}\right)}^2\hfill \\ \text{subject to}\hfill & 0<\frac{G_{\mathit{xy}}}{E_x}<1\hfill \\ \hfill & 0<\frac{G_{\mathit{xy}}}{E_y}<1\hfill \\ \hfill & L{B}_i\le x_i\le U{B}_i;i=1,2,3,\dots ,n\hfill \end{array}$$(6)

where f is the cost function of the design variable x = [x ₁, x ₂, x ₃, …, x _n ], ω _i is the evaluated natural frequencies using the Forward Fourier method, ${\overline{\mathit{\omega }}}_i$ is the benchmark or experimental natural frequencies, LB _i and UB _i are the lower and upper boundaries of the elastic properties, i.e. in-plane longitudinal elastic modulus (E _x ), in-plane transverse elastic modulus (E _y ), in-plane shear modulus (G _xy ) and the major Poisson’s ratio (v _xy ). As shown in Equation (6(6) $$\begin{array}{cc}\text{Minimise}\hfill & f\left(x\right)=\sum _{i=1}^N{\left(\frac{{\mathit{\omega }}_i-{\overline{\mathit{\omega }}}_i}{{\mathit{\omega }}_i}\right)}^2\hfill \\ \text{subject to}\hfill & 0<\frac{G_{\mathit{xy}}}{E_x}<1\hfill \\ \hfill & 0<\frac{G_{\mathit{xy}}}{E_y}<1\hfill \\ \hfill & L{B}_i\le x_i\le U{B}_i;i=1,2,3,\dots ,n\hfill \end{array}$$(6) ), two additional constraints are introduced to mitigate uncertainties as well as to avoid premature convergence.

4. Conventional algorithms

4.1. Genetic algorithm (GA)

GA was first developed by John Holland in the 1960s and since then, its potential advantages favour its widespread application in multiples fields. The basic idea of GA was inspired by the principle of genetic evolution in living organisms. Selection, crossover, mutation and elitism are the basic operators in genetic algorithms. An objective function is used to evaluate the fitness of each individual and the selection of individuals for further evolution or/and reproduction is made based on the fitness survival theory introduced by [31]. Crossover of genetic codes between pairs of chromosomes promotes local exploitation within a promising region while mutation of genetic codes in each chromosome encourages global exploration within a designated search area.

4.2. Ant colony optimization (ACO)

ACO was first invented by [32], emulating the foraging behaviour of real ants, in which they manage to build the shortest path between their settlement and the feeding source. Pheromone trails are used as a communication tool to construct the paths. Higher intensity of accumulated pheromones indicates higher possibility of prioritizing the path as the preferred one, and finally, generating the shortest solution path. The pheromone, τ _ij on the edge connecting cities i and j after iteration t is updated using equation as follows:(7) $${\mathit{\tau }}_{\mathit{ij}}\left(t+1\right)=\left(1-\mathit{\rho }\right).{\mathit{\tau }}_{\mathit{ij}}+\sum _{k=1}^m\mathrm{\Delta }{\mathit{\tau }}_{\mathit{ij}}^k$$(7) (8) $$\mathrm{\Delta }{\mathit{\tau }}_{\mathit{ij}}^k=\left\{\begin{array}{c}Q/L_k\hfill \\ 0\hfill \end{array}\right.\begin{array}{c}\text{if ant}k\text{uses edge}\left(i,j\right)\text{in its tour},\hfill \\ \text{otherwise,}\hfill \end{array}$$(8)

where ρ denotes the evaporation rate between iteration t and t + 1, m represents the number of ants, Q is a constant and L _k defines the length of tour built by kth ant.

The possibility of an ant k heading from city i to city j depends on the amount of accumulated pheromones, and hence, the probability function constituting the amount of pheromone is expressed as follows:(9) $$p_{\mathit{ij}}^k=\left\{\begin{array}{cc}\frac{{\mathit{\tau }}_{\mathit{ij}}^{\mathit{\alpha }}.{\mathit{\eta }}_{\mathit{ij}}^{\mathit{\beta }}}{{\sum }_{c_{\mathit{il}}\in N(s^p)}{\mathit{\tau }}_{\mathit{il}}^{\mathit{\alpha }}{\mathit{\eta }}_{\mathit{il}}^{\mathit{\beta }}}\hfill & \text{if}{c}_{\mathit{ij}}\in N\left(s^p\right)\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}\right.,$$(9)

where $N\left(s^p\right)$ is the set of cities available for kth ant, l denotes an unvisited city by kth ant, α and β represent the relative significance of the trail over the path information, ${\mathit{\eta }}_{\mathit{ij}}=1/d_{\mathit{ij}}$ with d _ij denoting the distance between cities i and j.

4.3. Particle swarm optimization (PSO)

PSO was first constructed by [33], mimicking the social character of animals (e.g. a flock of birds, a school of fish and a swarm of bees) to solve optimization problems. Unlike GA, PSO utilizes the information of local as well as the global best positions of particles to solve optimization problems. The formulations of PSO are defined as follows:(10) $${\overrightarrow{v}}_{k+1}=\overrightarrow{a}\otimes {\overrightarrow{v}}_k+{\overrightarrow{b}}_1\otimes {\overrightarrow{r}}_1\otimes \left({\overrightarrow{p}}_1-{\overrightarrow{x}}_k\right)+{\overrightarrow{b}}_2\otimes {\overrightarrow{r}}_2\otimes ({\overrightarrow{p}}_2-{\overrightarrow{x}}_k)$$(10) (11) $${\overrightarrow{x}}_{k+1}=\overrightarrow{c}\otimes {\overrightarrow{x}}_k+\overrightarrow{d}\otimes {\overrightarrow{v}}_{k+1}$$(11)

where is the updated velocity, ${\overrightarrow{v}}_{k+1}$ indicates element-by-element vector multiplication, k is the number of iteration, ${\overrightarrow{v}}_k$ is the previous velocity, $\overrightarrow{a}$ is the momentum factor, ${\overrightarrow{p}}_1$ is the previous best position, ${\overrightarrow{p}}_2$ is the global best position in the whole neighbourhood, ${\overrightarrow{b}}_1$ and ${\overrightarrow{b}}_2$ are the strength of attraction, ${\overrightarrow{x}}_{k+1}$ is the updated position, ${\overrightarrow{x}}_k$ is the previous position,$\overrightarrow{c}$ and $\overrightarrow{d}$ are the coefficients (normally unity). The assignment of input variables in the equations above decides the exploration and exploitation processes, depending on types of problems.

5. The proposed algorithm

The primary intention of combining GA, ACO and PSO is to minimize each other’s weaknesses and in the meantime, to promote respective strengths and therefore, transforming into a better hybrid algorithm. It should be highlighted that instead of eliminating the random initialization of GA, the effect of random initialization will be subdued by passing the products of GA through the ACO and PSO operators, whereby, in the composition of the proposed algorithm, ACO is assigned for exploratory search, while PSO is designated for exploitative search. The products of GA will be further processed by ACO and PSO to well distribute and organize the search coverage and thus, neutralize the biased effect of random initialization. Consequently, this will enhance the searching process. The collaboration between GA and ACO enhances the exploratory search capability such that ACO evades repeated tours from happening and the presence of pheromone trails probability function makes the exploration process to be more directive and organized. Meanwhile, the integration between GA and PSO improves the exploitative search capability as PSO is designated in only exploitative search to help GA since exploitation is comparatively the strength of PSO. In return, with the introduction of different operators, GA enhances the convergence rate and subdues the substantial dependency on the array of pheromones in ACO as well as further boosts the exploitative search in PSO. Two types of mutation operators are suggested in this algorithm, namely, two-point standard mutation and one-point refined mutation. In early iterations, standard mutation will be utilized collaboratively with the concept of unrepeated tour in ACO to evade local entrapment, while refined mutation will be used in later iterations to supplement the PSO exploitative search. The detailed procedures of the GA-ACO-PSO are described as follows.

5.1. Initialization

Number of generations, number of populations and number of dimensions (design variables) are specified accordingly. Before running the algorithm, input variables for ACO (e.g. pheromone initial intensity, number of paths, constant value) and input variables for PSO (e.g. momentum factor, strength of attraction) are specified appropriately. The lower and upper boundaries of each design variable are predefined. Once assignment of inputs is done, the algorithm is set to run. The algorithm begins with random generation of two populations, namely, the ant population and the particle population. It is noted that the number of ant population is equal to the number of particle population. These randomly generated solutions are then used to initialize the GA operators. Fitness value of each set of vector is then evaluated.

5.2. Selection, crossover, mutation, elitism

Selection is made based on the fitness values, depending on the types of problems (maximization or minimization). Vectors that possess better fitness values stand a higher possibility of getting selected for crossover operation. Selected vectors of ACO and PSO are then carried on to one-point inter-crossover. Thereafter, two-point mutation operation takes place for exploration purpose, followed by evaluation of fitness. In the stage of mutation, it consists of two different assignments: two-point standard mutation and one-point refined mutation. In early iterations, standard mutation is used for exploration reason. The concept of unvisited route in ACO is adopted to avoid repetition of tours, and thus, might lead to better convergence rate. Meanwhile, in later iterations when there is no unvisited route in the list or when the iteration exceeds half of its assigned value, standard mutation will be replaced by refined mutation. The objective of introducing refined mutation is to further improve the exploitation process which is mainly contributed by PSO. Before proceeding to variable updating, elitism is conducted from the set of vectors that comprises vectors of before and after evolution processes with the best individual remained at the top of the list and the worst at the bottom. The number of elite individuals that are remained for next iteration corresponds to the predefined number of populations. In the proposed algorithm, odd number elite individuals from the list with promising fitness values are eventually treated as the new populations of ACO. Comparison of the new population’s fitness with the previous one before undergoing evolution processes is made to ensure better exploration in the next iteration. If the new solution possesses better fitness, the previous solution will be replaced. Meanwhile, for the case of PSO, even number elite individuals from the elite group are used to update the particles best and global best positions, hence, initiating the algorithm in the next iteration. It is noted that the initial velocity vectors are set equal to zero and begin from the second iteration. After the replacement of PSO-evaluated solutions with the elite individuals, the velocity vectors are updated as well by the subtraction of the evaluated solutions from the elite solutions. This step should be taken into account for the reason of preserving the exactness of the velocity, corresponding to the distance between previous and current particles’ positions.

The formulations of crossover in the proposed algorithm are defined as follows:(12) $${\overrightarrow{x}}_{cross1}=[(rand1\times {\overrightarrow{x}}_1\left(1:co\right))+((1-rand1)\times {\overrightarrow{x}}_2\left(1:co\right)),(rand1\times {\overrightarrow{x}}_2\left(co+1:end\right))+((1-rand1)\times {\overrightarrow{x}}_1\left(co+1:end\right))]$$(12)

(13) $${\overrightarrow{x}}_{cross2}=[(rand2\times {\overrightarrow{x}}_2\left(1:co\right))+((1-rand2)\times {\overrightarrow{x}}_1\left(1:co\right)),(rand2\times {\overrightarrow{x}}_1\left(co+1:end\right))+((1-rand2)\times {\overrightarrow{x}}_2\left(co+1:end\right))]$$(13)

where ${\overrightarrow{x}}_{cross1}$ and ${\overrightarrow{x}}_{cross2}$ are the vectors after undergoing crossover process, rand1 and rand2 are random numbers ranging from 0 to 1, ${\overrightarrow{x}}_1$ and ${\overrightarrow{x}}_2$ denote the pair of vectors from ACO and PSO, and co indicates the crossover point.

Meanwhile, the equations of two-point standard mutation are expressed as follows:(14) $$x_{1new}=L{B}_1+\mathit{\alpha }\times (U{B}_1-L{B}_1)$$(14) (15) $$x_{2new}=L{B}_2+\mathit{\beta }\times (U{B}_2-L{B}_2)$$(15) (16) $$\mathit{\alpha }=(x_1-L{B}_1)/(U{B}_1-L{B}_1)$$(16)

(17) $$\mathit{\beta }=(x_2-L{B}_2)/(U{B}_2-L{B}_2)$$(17)

where x _1new and x _2new are the two mutated variables randomly selected from a vector $\overrightarrow{x}$, LB and UB are the predefined lower and upper boundaries of the variables, and x ₁ and x ₂ are the two variables before undergoing mutation process. In order to avoid any repeated tours in standard mutation, all the possible routes are constructed based on the permutation of different path number for each design variable. The route for each of the mutated design variables is determined and visited routes are deleted from the list. Repeated route will be replaced with another random pick of route from the list of unvisited routes. Furthermore, the concept of one-point refined mutation is defined as follows:

(18) $$\text{percent}=aa+rand\times (bb-aa)$$(18)

(19) $$lb=x\times \left(1-\text{percent}\right)$$(19)

(20) $$ub=x\times (1+\text{percent})$$(20)

(21) $$x_{\mathit{new}}=lb+rand\times (ub-lb)$$(21)

where aa and bb represent the lower and upper boundaries of percentage that need to be specified, percent indicates the randomly picked percentage within the specified range, x is one of the variables randomly selected from the best solution vector, lb and ub denote the lower and upper boundaries of variable, x and $x_{\text{new}}$ is the mutated variable. The concept of refined mutation is depicted in Figure 1, in which x ₂ has been picked for refined mutation. It is noted that the lower boundary or/and upper boundary in refined mutation will be set equal to the predefined boundaries in case that the predefined values are exceeded.

Figure 1. Concept of refined mutation.

5.3. Variable updating

The intensity of pheromone in ACO is updated using equations (7(7) $${\mathit{\tau }}_{\mathit{ij}}\left(t+1\right)=\left(1-\mathit{\rho }\right).{\mathit{\tau }}_{\mathit{ij}}+\sum _{k=1}^m\mathrm{\Delta }{\mathit{\tau }}_{\mathit{ij}}^k$$(7) -8(8) $$\mathrm{\Delta }{\mathit{\tau }}_{\mathit{ij}}^k=\left\{\begin{array}{c}Q/L_k\hfill \\ 0\hfill \end{array}\right.\begin{array}{c}\text{if ant}k\text{uses edge}\left(i,j\right)\text{in its tour},\hfill \\ \text{otherwise,}\hfill \end{array}$$(8) ), while for PSO, the velocities of particles are assigned randomly as well as the local and best positions of the particle are determined and inserted into equations (10(10) $${\overrightarrow{v}}_{k+1}=\overrightarrow{a}\otimes {\overrightarrow{v}}_k+{\overrightarrow{b}}_1\otimes {\overrightarrow{r}}_1\otimes \left({\overrightarrow{p}}_1-{\overrightarrow{x}}_k\right)+{\overrightarrow{b}}_2\otimes {\overrightarrow{r}}_2\otimes ({\overrightarrow{p}}_2-{\overrightarrow{x}}_k)$$(10) -11(11) $${\overrightarrow{x}}_{k+1}=\overrightarrow{c}\otimes {\overrightarrow{x}}_k+\overrightarrow{d}\otimes {\overrightarrow{v}}_{k+1}$$(11) ) to initiate the algorithm in the second iteration.

5.4. Solution construction

In the next iteration of ACO, the solutions are constructed using Equations (22(22) $$p_{\mathit{ij}}^k=\frac{{\mathit{\tau }}_{\mathit{ij}}}{{\sum }_{j=1}^n{\mathit{\tau }}_{\mathit{ij}}}$$(22) ) and (23(23) $$x_i^k=\frac{(pat{h}_{\mathit{number}}+rand-1)\times (U{B}_i-L{B}_i)}{\text{Number of paths}}+L{B}_i$$(23) ) as shown below. The concept of combined discrete and continuous probability distribution scheme is demonstrated in Figure 2 as shown below. The predefined range of the design variable is divided equally by the discrete total number of paths. Random pick of each design variable is then made within the region of respective selected paths. The ‘circles’ denote the ants or solutions with respective routes.

Figure 2. Concept of combined discrete and continuous probability distribution scheme.

Instead of using Equation (9(9) $$p_{\mathit{ij}}^k=\left\{\begin{array}{cc}\frac{{\mathit{\tau }}_{\mathit{ij}}^{\mathit{\alpha }}.{\mathit{\eta }}_{\mathit{ij}}^{\mathit{\beta }}}{{\sum }_{c_{\mathit{il}}\in N(s^p)}{\mathit{\tau }}_{\mathit{il}}^{\mathit{\alpha }}{\mathit{\eta }}_{\mathit{il}}^{\mathit{\beta }}}\hfill & \text{if}{c}_{\mathit{ij}}\in N\left(s^p\right)\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}\right.,$$(9) ) as mentioned above, solutions are probabilistically constructed based on the following probability Equation (22(22) $$p_{\mathit{ij}}^k=\frac{{\mathit{\tau }}_{\mathit{ij}}}{{\sum }_{j=1}^n{\mathit{\tau }}_{\mathit{ij}}}$$(22) ):

(22) $$p_{\mathit{ij}}^k=\frac{{\mathit{\tau }}_{\mathit{ij}}}{{\sum }_{j=1}^n{\mathit{\tau }}_{\mathit{ij}}}$$(22)

where τ _ij is the pheromone intensity for ith design variable at jth path. At initial stage, the values of pheromone intensity for each path are set equal to unity and thus, the probability of selecting either path is equal for the first iteration. The discrete path number is picked based on Equation (22(22) $$p_{\mathit{ij}}^k=\frac{{\mathit{\tau }}_{\mathit{ij}}}{{\sum }_{j=1}^n{\mathit{\tau }}_{\mathit{ij}}}$$(22) ) and the solution is randomly picked within a continuous range of a path number using Equation (23(23) $$x_i^k=\frac{(pat{h}_{\mathit{number}}+rand-1)\times (U{B}_i-L{B}_i)}{\text{Number of paths}}+L{B}_i$$(23) ) as shown below. It is noted that the use of this probability scheme is to mainly assist the exploration search in GA.

(23) $$x_i^k=\frac{(pat{h}_{\mathit{number}}+rand-1)\times (U{B}_i-L{B}_i)}{\text{Number of paths}}+L{B}_i$$(23)

where $x_i^k$ is the ith design variable of kth ant, path _number is the path number that is picked based on the probability of pheromone accumulation (initially, all the probability are set equal to unity), $\text{rand}$ denotes the random number ranges from 0 to 1, UB _i and LB _i are upper and lower boundaries of $i\text{th}$ design variable, and $\text{Number of paths}$ indicates the predefined total number of paths. Meanwhile, in PSO, the previous velocities and the previous solution vectors will be updated using Equations (10(10) $${\overrightarrow{v}}_{k+1}=\overrightarrow{a}\otimes {\overrightarrow{v}}_k+{\overrightarrow{b}}_1\otimes {\overrightarrow{r}}_1\otimes \left({\overrightarrow{p}}_1-{\overrightarrow{x}}_k\right)+{\overrightarrow{b}}_2\otimes {\overrightarrow{r}}_2\otimes ({\overrightarrow{p}}_2-{\overrightarrow{x}}_k)$$(10) ) and (11(11) $${\overrightarrow{x}}_{k+1}=\overrightarrow{c}\otimes {\overrightarrow{x}}_k+\overrightarrow{d}\otimes {\overrightarrow{v}}_{k+1}$$(11) ), thus, constituting new velocities and solution vectors in the next iteration. The proposed hybrid algorithm will continue operating until the maximum number of generations has taken place. For repeatability investigation, five independent runs are initiated. The complete course of the proposed hybrid algorithm is presented in Figure 3 and the pseudocode is shown in Figure 4.

Figure 3. The proposed hybrid GA-ACO-PSO algorithm.

Figure 4. Pseudocode of the proposed hybrid algorithm.

6. Results and discussions

6.1. Verification examples

The effectiveness of GA, PSO, ACO and the proposed hybrid algorithm in identifying the elastic properties of graphite epoxy plates and aluminium plates is investigated and compared. In the present study, Fourier method introduced by [28] is employed in the identification process as the Forward method. Three different sets of general boundary conditions are investigated, i. e. F-F-F-F, C-F-F-F and C-C-F-C (F: free; C: clamped). The accuracy of optimization algorithms is studied and compared in terms of minimized objective function values. The determined elastic properties are subsequently compared with the benchmark properties for quantitative verification. For repeatability study, five runs of simulations are executed and the standard deviations are evaluated. Furthermore, convergence study is conducted for different types of materials and different sets of boundary conditions.

The specifications of graphite epoxy plates and aluminium plates are taken from [34]. For graphite epoxy plates, the reference natural frequencies generated using Fourier method are retrieved from [28]. On the other hand, the reference natural frequencies of aluminium plates are generated using the mentioned published Fourier method with inputs of actual material properties taken from [34]. The dimensions and material properties of the plates are presented in Table 1. Meanwhile, the Fourier-generated reference frequencies for both the plates with different sets of boundary conditions are presented in Table 2.

Table 1. Dimensions and material properties of plates retrieved from [34].

Boundary condition	a, b (mm)	h (mm)	ρ (kgm^−3)	E x (GPa)	E y (GPa)	G xy (GPa)	v xy	v yx
Graphite epoxy plates
F-F-F-F		1.483
C-F-F-F	254.0	1.688	1584	127.9	10.27	7.312	0.22	0.0177
C-C-F-C		1.379
Aluminium plates
F-F-F-F
C-F-F-F	254.0	3.160	2770	72.4	72.4	28.0	0.33	0.33
C-C-F-C

Table 2. Reference natural frequencies of graphite epoxy plate and aluminium plate.

Boundary condition	f 1	f 2	f 3	f 4	f 5	f 6
Graphite epoxy plates
F-F-F-F	51.76	60.18	121.4	165.9	212.6	228.6
C-F-F-F	38.04	52.58	105.3	217.2	238.4	257.4
C-C-F-C	70.96	167.5	219.7	298.1	314.8	432.2
Aluminium plates
F-F-F-F	163.02	230.03	299.34	419.94	419.94	741.04
C-F-F-F	42.03	102.98	255.05	331.39	374.45	658.74
C-C-F-C	290.11	489.01	767.24	938.59	984.95	1431.31

The identification process involves minimization of the objective function and the input settings for GA, ACO, PSO and hybrid algorithms are listed in Table 3. The number of evaluations is set constant for all the involved algorithms by fixing the number of generations and number of populations for the reason of equal and fair comparison. In the present study, the number of evaluations includes those from initialization stage as well. Basically, the number of evaluations can be interpreted as the number of times that the minimization of the objective function has taken place. The number of evaluations for GA, ACO and PSO can be calculated based on Equation (24(24) $$\text{Number of evaluations}=\text{Number of generations}\left(\text{Number of populations}\right)+\text{Number of initial populations}$$(24) ), while for the proposed hybrid algorithm, the number of evaluations can be determined using Equation (25(25) $$\text{Number of evaluations}=2(\text{Number of generations}\left(\text{Number of populations}\right)+\text{Number of initial populations})$$(25) ).(24) $$\text{Number of evaluations}=\text{Number of generations}\left(\text{Number of populations}\right)+\text{Number of initial populations}$$(24) (25) $$\text{Number of evaluations}=2(\text{Number of generations}\left(\text{Number of populations}\right)+\text{Number of initial populations})$$(25)

Table 3. Input settings for GA, ACO, PSO and hybrid algorithm.

	Input variables	Setting 1
		GA	ACO	PSO	Hybrid
1.	Number of generations	20	20	20	20
2.	Number of populations	40	40	40	20
3.	Number of evaluations	840	840	840	840
4.	Number of dimensions (design variables)	4	4	4	4
5.	ACO populations	–	40	–	10
6.	Initial intensity of pheromone for every vectors	–	1	–	1
7.	Number of paths	–	30	–	25
8.	Evaporation rate	–	0	–	0
9.	Constant	–	10,000	–	10,000
10.	PSO populations	–	–	40	10
11.	Strength of attraction, ${\overrightarrow{b}}_1,{\overrightarrow{b}}_2$	–	–	1.5, 1.5	1.05, 1.05
12.	Momentum factor, $\overrightarrow{a}$	–	–	[0.9,0.4]; linspace = 20	[0.5,0.1]; linspace = 20
13.	Number of crossover	20	–	–	10
14.	Number of standard (or refined) mutation	20	–	–	10 (10)
15.	Initialization of refined mutation after kth generation	–	–	–	10th
16.	The percentage range of design variables in refined mutation, aa & bb	–	–	–	aa = 0.02; bb = 0.05;

In the present study, four dimensions are involved, e.g. in-plane longitudinal elastic modulus (E _x ), in-plane transverse elastic modulus (E _y ), in-plane shear modulus (G _xy ) and the major Poisson’s ratio (v _xy ). The optimal number of path used in ACO and hybrid algorithm that yield the best results in material identification is found to be 30 and 25, respectively. It is noted that higher number of path leads to higher degree of exploitative search. Hence, it can be seen that in hybrid algorithm, the contribution of ACO operator is mainly focalized on the exploratory search. On the other hand, the constant is a fixed number involved in pheromone intensity updating in ACO. It is specified based on the type of objective function and the optimal constant used in this context is 10,000 for both ACO and hybrid algorithm. Besides, the strength of attraction and the momentum factor in PSO are specified referring to published article; while for hybrid algorithm, those parameters are set lower so that the PSO operator in the algorithm could concentrate more on the exploitative search. The command linspace in PSO, is used to generates a linearly spaced vector within a specified range in descending order, for instance, when linspace = 20, from 0.9 to 0.4, there will be 20 linearly spaced interval within the specified range. In this study, the linspace is specified corresponding to the number of generations. Moreover, number of crossover is the number of selected individual to undergo crossover process and similarly, number of mutation is the number of selected individual to undergo mutation process.

6.1.1. Accuracy

The indicators evaluating the effectiveness of these algorithms include the value of minimization (accuracy), the standard deviation (repeatability), the rate of convergence and computational time (convergence) as well as the robustness. Figures 5 and 6 demonstrate the comparison of minimized objective function values achieved by GA, ACO, PSO and the proposed hybrid algorithm for both graphite epoxy and aluminium plates, respectively. It is noted that for visual qualitative comparison reason, the graphs shown in Figures 5 and 6 have been enlarged and the maximum y-axis limit does not represent the actual maximum values of those algorithms. From Figure 5, it can be seen that PSO, GA and the proposed hybrid algorithm yield considerably better outcomes compared to ACO. The results suggest that PSO and GA can be the probable alternatives as it can be observed that the values are comparatively close but marginally inferior to those of the proposed hybrid algorithm. In comparison of the four algorithms, the proposed hybrid algorithm yields the best results for every set of boundary conditions with values of 2.8081E−08, 1.2745E−07 and 3.5636E−08, as shown in Tables 4–6, respectively. This consequently verifies the accuracy of the proposed algorithm. Furthermore, the absolute percentage errors of the determined elastic properties with respect to benchmark properties using the proposed hybrid algorithm are found to be below 3% and this certainly acknowledges the reliability of the algorithm. From Figure 6, it can be observed that ACO performs the worst, followed by GA and PSO, while hybrid algorithm demonstrates relatively the best results. Still, in the study of aluminium plate, the proposed hybrid algorithm appears to be the relatively the best choice with minimized values of 1.5579E−08, 2.4244E−08 and 7.8573E−08 for three different sets of boundary conditions, as presented in Tables 7–9, respectively. The improvement percentage of the identified elastic properties when using the proposed algorithm with respect to those using conventional algorithms is presented in Table 10. It shows that the accuracy of the determined E _x and E _y using the proposed hybrid algorithm in comparisons to those of using other algorithms is improved by 0.3836% and 2.2911% at most for graphite epoxy and aluminium plates, respectively. Furthermore, the evaluated in-plane shear modulus of graphite epoxy and aluminium plates utilizing the proposed hybrid algorithm is found to be more accurate by at most 0.8958% and 1.5589%, respectively, when comparing to those of using traditional algorithms. Significant improvements can as well be seen in the identified Poisson’s ratio of graphite epoxy plate (by 19.410% at most) and aluminium plate (by 7.5151% at most) when applying the proposed hybrid algorithm. In an overall view, the proposed hybrid algorithm is proven to be slightly more superior to the others in terms of accuracy.

Figure 5. Comparison of minimized value for graphite epoxy plate using different algorithms.

Figure 6. Comparison of minimized value for aluminium plate using different algorithms.

Table 4. Elastic constants obtained for graphite epoxy plate with F-F-F-F boundary condition.

	E x (GPa)	E y (GPa)	G xy (GPa)	v xy	Obj. Fun.
PSO	127.8802 [0.0155]	10.2697 [0.0029]	7.3092 [0.0383]	0.2067 [6.0455]	2.5121E−07
GA	127.8441 [0.0437]	10.2668 [0.0312]	7.3117 [0.0041]	0.2316 [5.2727]	2.1230E−07
ACO	128.4253 [0.4107]	10.2508 [0.0182]	7.3072 [0.0656]	0.1756 [20.182]	6.9808E−06
Hybrid	127.8654 [0.0271]	10.2696 [0.0039]	7.3134 [0.0191]	0.2183 [0.7727]	2.8081E−08
Benchmark	127.9	10.27	7.312	0.22	–

Note: [n]Absolute percentage error.

Table 5. Elastic constants obtained for graphite epoxy plate with C-F-F-F boundary condition.

	E x (GPa)	E y (GPa)	G xy (GPa)	v xy	Obj. Fun.
PSO	127.8159 [0.0658]	10.2619 [0.0789]	7.3263 [0.1956]	0.2368 [7.6364]	3.2977E−07
GA	127.8445 [0.0434]	10.2391 [0.3009]	7.3410 [0.3966]	0.1957 [11.045]	1.3924E−06
ACO	128.2228 [0.2524]	10.3010 [0.3019]	7.2463 [0.8985]	0.2630 [19.545]	8.0426E−06
Hybrid	127.9189 [0.0148]	10.2781 [0.0789]	7.3118 [0.0027]	0.2140 [2.7273]	1.2745E−07
Benchmark	127.9	10.27	7.312	0.22	–

Note: [n]Absolute percentage error.

Table 6. Elastic constants obtained for graphite epoxy plate with C-C-F-C boundary condition.

	E x (GPa)	E y (GPa)	G xy (GPa)	v xy	Obj. Fun.
PSO	127.6372 [0.2055]	10.2642 [0.0565]	7.2971 [0.2038]	0.2362 [7.3636]	7.9172E−07
GA	128.0523 [0.1191]	10.2709 [0.0088]	7.2975 [0.1983]	0.2186 [0.6364]	1.3538E−07
ACO	128.3453 [0.3482]	10.2966 [0.2590]	7.2949 [0.2339]	0.2159 [1.8636]	5.2515E−06
Hybrid	127.8287 [0.0557]	10.2688 [0.0117]	7.3164 [0.0602]	0.2203 [0.1364]	3.5636E−08
Benchmark	127.9	10.27	7.312	0.22	–

Note: [n]Absolute percentage error.

Table 7. Elastic constants obtained for aluminium plate with F-F-F-F boundary condition.

	E x (GPa)	E y (GPa)	G xy (GPa)	v xy	Obj. Fun.
PSO	72.4631 [0.0872]	72.3508 [0.0680]	28.0034 [0.0121]	0.3330 [0.9091]	3.0826E−06
GA	72.4607 [0.0838]	72.3786 [0.0296]	27.9958 [0.0150]	0.3301 [0.0303]	1.1039E−07
ACO	72.7990 [0.5511]	72.1331 [0.3686]	27.9286 [0.2550]	0.3344 [1.3333]	8.7278E−06
Hybrid	72.4002 [0.0003]	72.4211 [0.0291]	28.0019 [0.0068]	0.3299 [0.0303]	1.5579E−08
Benchmark	72.4	72.4	28	0.33	–

Note: [n]Absolute percentage error.

Table 8. Elastic constants obtained for aluminium plate with C-F-F-F boundary condition.

	E x (GPa)	E y (GPa)	G xy (GPa)	v xy	Obj. Fun.
PSO	72.3556 [0.0613]	73.0242 [0.8622]	28.0288 [0.1029]	0.3158 [4.3030]	4.4058E−06
GA	72.4287 [0.0396]	72.8125 [0.5698]	28.0457 [0.1632]	0.3159 [4.2727]	4.8509e−06
ACO	72.4018 [0.0025]	70.7281 [2.3091]	27.9373 [0.2239]	0.3553 [7.6667]	2.5178E−05
Hybrid	72.4057 [0.0079]	72.4130 [0.0180]	27.9956 [0.0157]	0.3295 [0.1515]	2.4244e−08
Benchmark	72.4	72.4	28	0.33	–

Note: [n]Absolute percentage error.

Table 9. Elastic constants obtained for aluminium plate with C-C-F-C boundary condition.

	E x (GPa)	E y (GPa)	G xy (GPa)	v xy	Obj. Fun.
PSO	72.1016 [0.4122]	72.2972 [0.1420]	27.7662 [0.8350]	0.3390 [2.7273]	3.6406e−07
GA	72.2092 [0.2635]	72.2842 [0.1599]	27.8983 [0.3632]	0.3355 [1.6667]	4.0050e−07
ACO	73.5411 [1.5761]	72.7066 [0.4235]	28.6222 [2.2221]	0.3028 [8.2424]	4.3167e−06
Hybrid	72.3456 [0.0751]	72.3924 [0.0105]	27.8143 [0.6632]	0.3342 [1.2727]	7.8573e−08
Benchmark	72.4	72.4	28	0.33	–

Note: [n]Absolute percentage error.

Table 10. Improvement percentage of identified elastic properties using the proposed algorithm with respect to those of using other algorithms.

	E x (%)	E y (%)	G xy (%)	v xy (%)
Graphite epoxy	0.0166–0.3836	0.0143–0.2473	0.0192–0.8958	0.5000–19.410
Aluminium	0.0317–1.501	0.0005–2.2911	0.0053–1.5589	0.3940–7.5152

6.1.2. Repeatability

Repeatability can be defined as the consistency of the outcomes yielded by an algorithm. It is evaluated in terms of standard deviation, in which a higher value of standard deviation indicates lower repeatability. From Figures 7 and 8, the proposed hybrid algorithm and GA are found to be the better options as compared to PSO and ACO. As a whole, the proposed hybrid algorithm performs the best with standard deviations of 7.7340E−08, 1.0473E−07, 2.3007E−07 for graphite epoxy plate with three different sets of boundary conditions, as demonstrated in Tables 11–13, respectively, as well as the standard deviations for aluminium plate with three different sets of boundary conditions are found to be 2.7751E−07, 1.9999E−07, 1.8245E−07, as presented in Tables 14–16, respectively. From Table 17, the proposed hybrid algorithm can be seen leading the rest of the algorithms in terms of repeatability by at most 6.5410E−04 and 8.7205E−05 in difference of standard deviation for graphite epoxy and aluminium plates, respectively.

Figure 7. Comparison of repeatability for graphite plate using different algorithms.

Figure 8. Comparison of repeatability for aluminium plate using different algorithms.

Table 11. Repeatability study using different type of algorithms for graphite epoxy plate with F-F-F-F boundary condition.

	Best	Average	Worst	Std. dev.
PSO	2.5121E−07	3.3108E−04	1.6394E−03	6.5418E−04
GA	2.1230E−07	1.8251E−06	5.6053E−06	2.0889E−06
ACO	6.9808E−06	4.7169E−04	1.5170E−03	5.5143E−04
Hybrid	2.8081E−08	1.6298E−07	2.3106E−07	7.7340E−08

Table 12. Repeatability study using different type of algorithms for graphite epoxy plate with C-F-F-F boundary condition.

	Best	Average	Worst	Std. dev.
PSO	3.2977E−07	1.0014E−06	1.9707E−06	5.4871E−07
GA	1.3924E−06	5.2336E−06	1.2572E−05	4.1214E−06
ACO	8.0426E−06	2.3155E−04	5.4768E−04	1.8497E−04
Hybrid	1.2745E−07	2.2914E−07	3.8962E−07	1.0473E−07

Table 13. Repeatability study using different type of algorithms for graphite epoxy plate with C-C-F-C boundary condition.

	Best	Average	Worst	Std. dev.
PSO	7.9172E−07	3.4247E−05	1.5742E−04	6.1708E−05
GA	1.3538E−07	3.1563E−06	1.3302E−05	5.0936E−06
ACO	5.2515E−06	4.9431E−04	1.6992E−03	6.2138E−04
Hybrid	3.5636E−08	2.4207E−07	6.7863E−07	2.3007E−07

Table 14. Repeatability study using different type of algorithms for aluminium plate with F-F-F-F boundary condition.

	Best	Average	Worst	Std. dev.
PSO	3.0826E−06	1.8734E−05	5.2624E−05	1.8059E−05
GA	1.1039E−07	7.5095E−06	2.6547E−05	9.7445E−06
ACO	8.7278E−06	1.3182E−04	2.5658E−04	8.7483E−05
Hybrid	1.5579E−08	4.0034E−07	8.7327E−07	2.7751E−07

Table 15. Repeatability study using different type of algorithms for aluminium plate with C-F-F-F boundary condition.

	Best	Average	Worst	Std. dev.
PSO	4.4058E−06	2.2551E−05	8.1018E−05	2.9593E−05
GA	4.8509E−06	1.3303E−05	2.2814E−05	7.8665E−06
ACO	2.5178E−05	9.2809E−05	2.4608E−04	7.9389E−05
Hybrid	2.4244E−08	2.2518E−07	5.2392E−07	1.9999E−07

Table 16. Repeatability study using different type of algorithms for aluminium plate with C-C-F-C boundary condition.

	Best	Average	Worst	Std. dev.
PSO	3.6406E−07	1.3024E−05	4.3623E−05	1.6081E−05
GA	4.0050E−07	6.6893E−07	8.8026E−07	2.0707E−07
ACO	4.3167E−06	2.4678E−05	8.4013E−05	2.9836E−05
Hybrid	7.8573E−08	2.6036E−07	5.9360E−07	1.8245E−07

Table 17. Difference of standard deviation of the proposed hybrid algorithm with respect to those of other algorithms.

	F-F-F-F	C-F-F-F	C-C-F-C
Graphite epoxy	2.0116E−04 to 6.5410E−04	4.4398E−07 to 1.8487E−04	4.8635E−06 to 6.2115E−04
Aluminium	9.4670E−06 to 8.7205E−05	7.6665E−06 to 7.9189E−05	2.4620E−08 to 2.9654E−05

6.1.3. Convergence rate and computational time

Convergence rate is one of the key factors determining the reliability of an algorithm. It is usually regarded as the higher the rate of convergence, the better the algorithm, provided that the outcomes converge and fall within the acceptable range. In the present section, the manner of how fast the algorithm converges and attains its lowest minimized value in sequence of number of iterations is investigated and discussed. Figures 9–11 demonstrate the convergence graphs of the four algorithms for graphite epoxy plate with F-F-F-F, C-F-F-F and C-C-F-C boundary conditions, respectively, while, Figures 12–14 depict the convergence behaviours of the four algorithms for aluminium plate with F-F-F-F, C-F-F-F and C-C-F-C boundary conditions, respectively. Based on Figure 9, it can be observed that from iteration 1 to 3, there is a drastic reduction in the objective function value of the proposed hybrid algorithm before it is further reduced in iteration 6. In view of other algorithms, the initial values are comparatively smaller and the reduction of objective function values is not as much as that of the proposed hybrid algorithm. Meanwhile, for C-F-F-F boundary condition, the proposed hybrid algorithm exhibits quite similar rate of convergence as those of GA and ACO starting from iteration 3, as depicted in Figure 10. The proposed hybrid algorithm experiences a great drop in the objective function value from iteration 1 to 3. As shown in Figure 11, it can be observed that the proposed hybrid algorithm performs the best, in which it yields the least minimized objective function value at the highest rate among all the algorithms. Apart from the study of graphite epoxy, Figure 12 shows the convergence plot for aluminium plate with F-F-F-F boundary condition. As expected, the hybrid algorithm exhibits as good rate of convergence as that of ACO, while those of GA and PSO come after. In Figure 13, the hybrid algorithm once again demonstrates its great dominance in rate of convergence over the other algorithms, subjecting to a great decrease in the objective function values after the first iteration. For C-C-F-C boundary condition shown in Figure 14, the hybrid algorithm surpasses the other algorithms in terms of convergence rate in early iterations, however in iteration 6, PSO converges slightly sooner than that of the hybrid algorithm. Despite, the hybrid algorithm yields considerably good results. In overall, it can be observed that in early iterations, the proposed hybrid algorithm converges to its solution at a comparatively higher rate. This phenomenon occurs due to the equal contribution of work load from both the ACO and PSO operators, in which ACO operator is assigned for exploratory search, while PSO operator is designated for exploitative search. The biased effect due to random initialization of GA operator is compromised via the use of ACO and PSO operators and in the next iteration, those processed outcomes, whereby ones from the products of exploration search (ACO) and the other from the products of exploitative search (PSO) are subsequently merged and evolved again via the use of GA operators. Therefore, at the end, the rate of minimization can be seen on the rise in comparison to those of other involved algorithms due to well distribution of search coverage.

Figure 9. Convergence comparison of graphite epoxy plate for boundary condition F-F-F-F.

Figure 10. Convergence comparison of graphite epoxy plate for boundary condition C-F-F-F.

Figure 11. Convergence comparison of graphite epoxy plate for boundary condition C-C-F-C.

Figure 12. Convergence comparison of aluminium plate for boundary condition F-F-F-F.

Figure 13. Convergence comparison of aluminium plate for boundary condition C-F-F-F.

Figure 14. Convergence comparison of aluminium plate for boundary condition C-C-F-C.

Technically, algorithm that consumes much computational time can be referred to as computationally expensive (high computational cost). There is a slight difference between convergence rate and computational time in the present study, in which, convergence rate refers to the ability of the algorithm to converge involving lower number of evaluations, whereby, computational time refers to the total time consumed by the algorithm to converge to its solution. In this context, the number of evaluations for each algorithm is fixed constant for equal comparison purposes and the computational time of each algorithm for graphite epoxy and aluminium plates are presented in Figures 15 and 16, respectively. It is noted that the computational time is dependent on the specifications of the computing device, in which, in the present study, a desktop with Intel(R) Core (TM) i5–2400 CPU @ 3.10 GHz processor is used to study the computational time of the all the involved algorithm for two different plates with F-F-F-F boundary condition. The results show that under a fixed number of evaluations, the computational time of the involved algorithms is almost similar with time difference not exceeding 12 s. The proposed hybrid algorithm is found consuming the most computational time (3872s) if compared to GA (3860s), ACO (3867s) and PSO (3864s) for graphite epoxy plate with F-F-F-F boundary condition. In comparison to the least time-consuming GA, the proposed algorithm is found to have lagged behind by 12 s at most, which is equivalent to 0.3109% in difference. Despite, the shortcoming of the proposed algorithm in this aspect is in fact offset by its relatively large accuracy improvement, as shown in Table 10, garnering 0.3836, 0.2473, 0.8958 and 19.41% of improvement at most for E _x , E _y , G _xy and v _xy , respectively. Furthermore, in the study of aluminium plate with F-F-F-F boundary condition, the proposed hybrid algorithm is also found tailing the rest of the algorithm in terms of computational time (3928s), in which, GA, ACO and PSO consumes 3922s, 3926s and 3924s, respectively. The time difference between the proposed hybrid algorithm and the least time-consuming GA is found to be around 6 s, indicating 0.1530% in difference with respect to the time consumed by GA. Despite, comparing to the accuracy improvement of the proposed hybrid algorithm, the drawback of the proposed hybrid algorithm in the aspect of computational time is as well overshadowed by its relatively more significant accuracy improvement, where, the accuracy of the identified E _x , E _y , G _xy and v _xy using the proposed hybrid algorithm are improved by at most 1.501, 2.2911, 1.5589 and 7.5152%, respectively, with respect to those of using conventional algorithms, as demonstrated in Table 10. In a general prospect, relatively rapid convergence of the hybrid algorithm in sequence of iterations, its high repeatability as well as its excellent capability in searching for global minimum are said to have compensated its shortcoming in the aspect of computational time.

Figure 15. Comparison in computational time of graphite epoxy plate.

Figure 16. Comparison in computational time of aluminium plate.

6.1.4. Robustness

In real cases, experimental data often comes with uncertainties and these uncertainties are predominantly due to the presence of various errors, i.e. prediction errors, measurement errors and implementation errors. In the present study, implementation and measurement errors are the emphasis since these errors might incur during the acquisition of vibration data. Uncertainties can be reduced experimentally and handled numerically. To minimize uncertainties incurred during experimental measurement, several important issues pertaining to experimental set-up should be taken into serious attention before conducting the experiment. These include the identification of the boundary condition of the real structure, the suitability of impact location and impact grid, the hardness suitability of impact hammer tip as well as the suitability of block size and sampling rate. By addressing these issues appropriately, the level of uncertainties incurred during experiment could greatly be reduced. It is noteworthy that these issues are in fact dependent on the specimen dimensions as well as the resonant frequency range of the specimen. Even though several precautionary steps have been strictly followed during experiment, it is learnt that in reality, the presence of uncertainties in the measured data are inevitable; therefore, it is necessary to curb these uncertainties numerically. In numerical evaluation, in order to accommodate these uncertainties as well as to avoid premature convergence, constraints should be introduced, as defined in Equation (6(6) $$\begin{array}{cc}\text{Minimise}\hfill & f\left(x\right)=\sum _{i=1}^N{\left(\frac{{\mathit{\omega }}_i-{\overline{\mathit{\omega }}}_i}{{\mathit{\omega }}_i}\right)}^2\hfill \\ \text{subject to}\hfill & 0<\frac{G_{\mathit{xy}}}{E_x}<1\hfill \\ \hfill & 0<\frac{G_{\mathit{xy}}}{E_y}<1\hfill \\ \hfill & L{B}_i\le x_i\le U{B}_i;i=1,2,3,\dots ,n\hfill \end{array}$$(6) ). In the present study, two additional inequality constraints are defined to mitigate the presence of uncertainties, in which, the ratio of in-plane shear modulus (G _xy ) to in-plane longitudinal elastic modulus (E _x ) as well as the ratio of in-plane shear modulus (G _xy ) to in-plane transverse elastic modulus (E _y ) are restricted to be within the range of 0 to 1. Furthermore, since the range of elastic properties is unknown in practical applications, the search boundaries of elastic properties can initially be set large enough to encompass possible search areas. In addressing uncertainties in the measured results, the influence of these errors/uncertainties on the end outcomes are studied and presented. Three different levels of errors, i.e. 1, 5 and 10% are artificially introduced into the ‘measured’ natural frequencies and these ‘erroneous’ natural frequencies are subsequently used as the reference parameters in the proposed hybrid algorithm. The reason of introducing three levels of errors is to investigate the influence of increasing errors on the absolute percentage difference of the identified elastic properties. The error level is set up to 10% considering it is in the worst-case scenario. Apart, this investigation is extended to graphite epoxy plate and aluminium plate to study the effects of the types of plate on the absolute percentage difference of identified elastic properties in the presence of errors. It is noted that the absolute percentage difference is evaluated on the basis of the absolute percentage difference between the identified elastic properties and the actual reference elastic properties. From Tables 18 and 19, the identified elastic moduli (E _x and E _y ), as well as the identified in-plane shear modulus (G _xy ) experience elevating deviations as the error level rises from 1% to 10% with the absolute percentage difference ranging from 0.01% to 8%. On the other hand, the evaluated Poisson’s ratio (v _xy ) fluctuates within the range of 0.02% to 13% in the presence of increasing level of error. The occurrence of this phenomenon could possibly due to the insignificant influence of natural frequency on the identified Poisson’s ratio, v _xy . During initialization of numerical evaluation, the value of the Poisson’s ratio is first randomly assigned. Because of the insignificant influence of natural frequency on the Poisson’s ratio, great changes in natural frequency induce little changes in the Poisson’s ratio, causing the identified Poisson’s to be locally trapped near its initial value. Therefore, it can be observed that the identified Poisson’s ratio fluctuates and deviates more from its actual value if compared to the other elastic constants. In overall, the proposed hybrid algorithm is proven to be able to produce reliable identification results in the presence of errors, regardless of types of material as well as boundary conditions.

Table 18. Identified elastic properties of graphite epoxy plate in the presence of increasing errors.

Reference	E x (GPa)	E y (GPa)	G xy (GPa)	v xy
	127.9	10.27	7.312	0.22
F-F-F-F
1%	128.1903 [0.2270]	10.2659 [0.0399]	7.3114 [0.0082]	0.2470 [12.2727]
5%	129.1287 [0.9607]	10.1255 [1.4070]	7.3225 [0.1436]	0.2209 [0.4091]
10%	129.7710 [1.4629]	10.7189 [4.3710]	7.2655 [0.6359]	0.2443 [11.0455]
C-F-F-F
1%	127.2377 [0.5178]	10.3559 [0.8364]	7.2500 [0.8479]	0.2267 [3.0455]
5%	126.4558 [1.1292]	10.5217 [2.4508]	7.1893 [1.6781]	0.2120 [3.6364]
10%	124.2640 [2.8428]	10.9467 [6.5891]	7.5320 [3.0088]	0.2205 [0.2273]
C-C-F-C
1%	127.4838 [0.3254]	10.1851 [0.8267]	7.3196 [0.1039]	0.2282 [3.7273]
5%	123.1655 [3.7017]	10.4110 [1.3729]	7.0117 [4.1069]	0.2293 [4.2273]
10%	118.9493 [6.9982]	10.4508 [1.7605]	7.6814 [5.0520]	0.2344 [6.5455]

Note: [n]Absolute percentage error.

Table 19. Identified elastic properties of aluminium plate in the presence of increasing errors.

Reference	E x (GPa)	E y (GPa)	G xy (GPa)	v xy
	72.4	72.4	28	0.33
F-F-F-F
1%	72.5913 [0.2642]	72.4991 [0.1369]	28.0085 [0.0304]	0.3290 [0.3030]
5%	71.1711 [1.6974]	72.2650 [0.1865]	28.2343 [0.8368]	0.3190 [3.3333]
10%	73.9361 [2.1217]	77.4832 [7.0210]	28.2909 [1.0389]	0.2947 [10.6970]
C-F-F-F
1%	72.4563 [0.0778]	71.9042 [0.6848]	27.7964 [0.7271]	0.3524 [6.7879]
5%	71.1265 [1.7590]	71.5846 [1.1262]	27.6695 [1.1804]	0.3626 [9.8788]
10%	70.3329 [2.8551]	68.9921 [4.7070]	28.5882 [2.1007]	0.2999 [9.1212]
C-C-F-C
1%	72.7744 [0.5171]	72.4439 [0.0606]	27.9578 [0.1507]	0.3243 [1.7273]
5%	71.0602 [1.8506]	71.9986 [0.5544]	29.7598 [6.2850]	0.3212 [2.6667]
10%	69.8662 [3.4997]	67.1080 [7.3094]	30.2124 [7.9014]	0.3091 [6.3333]

Note: [n]Absolute percentage error.

6.2. Application examples

Unidirectional laminated carbon epoxy composite plate, unidirectional laminated glass epoxy composite plate and aluminium composite panel (ACP) are utilized to investigate the applicability of the suggested algorithm. The dimensions and the experimental resonant frequencies of the plates are shown in Tables 20 and 21, respectively. From Table 21, it can be seen that eight, seven and five modes of fundamental natural frequencies are involved in the identification of elastic properties of carbon epoxy composite plate, glass epoxy composite plate and ACP, respectively.

Table 20. Dimensions of carbon epoxy composite plate, glass epoxy composite plate and aluminium composite panel (ACP).

Composite plates	Length (mm)	Width (mm)	Thickness (mm)	Density (kg m^−3)
Carbon epoxy	207.5	207.5	2.00	1535
E-glass epoxy	268.0	150.0	1.43	1932
ACP	300.0	300.0	3.00	1520

Table 21. Measured natural frequencies of carbon epoxy composite plate, glass epoxy composite plate and aluminium composite panel (ACP) with F-F-F-F boundary condition.

Composite plates	f 1	f 2	f 3	f 4	f 5	f 6	f 7	f 8
Carbon epoxy [35]	97.0	123.8	237.4	341.0	458.0	502.6	541.2	653.0
Glass epoxy [37]	72.3	97.6	176.7	197.4	238.3	277.5	348.0	–
ACP	52.7	82.1	107.0	142.0	149.0	–	–	–

The details of carbon epoxy composite plate are taken from [35], while the details of glass epoxy are taken from [36,37]. For ACP, EMA is conducted on real specimen to experimentally determine its natural frequencies, which to be used in the identification process. The experimental set-up is demonstrated in Figures 17 and 18. As seen from the figures, an impact hammer is used to strike the plate, and the response is received by a mounted accelerometer which is connected to a data acquisition system. The position of the accelerometer is fixed at a particular point and the impact hammer is roved from point 1 to point 25. Five knocks are imparted on each point and the average response is recorded. The collected data are then processed to obtain the plate’s experimental natural frequencies which to be used in the identification process.

Figure 17. Experimental set-up.

Figure 18. Predefined nodes and grids on ACP.

Only the proposed hybrid algorithm is adopted into application since it has been proven to be the best choice among all the four algorithms in the aspects of accuracy, repeatability and convergence in previous verification section. Furthermore, only F-F-F-F boundary condition is studied in these examples. It is notable that the closeness of experimental and evaluated natural frequencies signifies the accuracy as well as the validity of the determined elastic properties. Hence, it can be interpreted as smaller error between the experimental and evaluated natural frequencies leads to more accurate identified elastic properties. The input settings for the hybrid algorithm used in application examples are similar to those of verification examples as listed in Table 3.

The material properties of carbon epoxy composite plate, glass epoxy composite plates and aluminium composite panel (ACP) evaluated using the proposed hybrid algorithm are presented in Tables 22–24, respectively. In the case of carbon epoxy composite plate, the minimized objective function values are found to be consistent throughout the five runs of simulations with the least minimized value of 6.7995E−05. The best evaluated in-plane elastic moduli (E _x and E _y ), in-plane shear modulus (G _xy ) and the major Poisson’s ratio (v _xy ) are discovered to be 170.5685 GPa, 10.3249 GPa, 6.0626 GPa and 0.3047, respectively. In the case of glass epoxy composite plate, the best values of the evaluated material properties are 44.2024 GPa, 16.9951 GPa, 7.2141 GPa and 0.3200 for the in-plane elastic moduli (E _x and E _y ), in-plane shear modulus (G _xy ) and the major Poisson’s ratio (v _xy ), respectively, with the least minimized value of 2.0609E−05. In the application of aluminium composite panel (ACP), the best identified in-plane elastic moduli (E _x and E _y ), in-plane shear modulus (G _xy ), and the major Poisson’s ratio (v _xy ) are 12.6805 GPa, 10.0935 GPa, 3.4892 GPa and 0.3352, respectively, yielding the least minimized value of 9.1873E−04. From Table 25, it can be observed that in overall, the absolute percentage errors of the evaluated natural frequencies with respect to the experimental natural frequencies are less than 2%, where the maximum absolute percentage error for carbon epoxy is 1.2410%, while for glass epoxy, the maximum absolute percentage error is 0.2939% as well as the maximum absolute percentage error for ACP is 1.6851%.

Table 22. The identified material properties of carbon epoxy composite plate.

	Run 1	Run 2	Run 3	Run 4	Run 5	Average
E x (GPa)	169.6889	170.0624	170.5685	169.7381	170.6120	170.1340
E y (GPa)	10.3860	10.3592	10.3249	10.3453	10.3675	10.3566
G xy (GPa)	6.1265	6.0641	6.0626	6.0721	6.1566	6.0964
v xy	0.3073	0.3409	0.3047	0.3356	0.3246	0.3226
Obj. Fun.	2.1457E−04	7.5633E−05	6.7995E−05	8.2959E−05	2.2989E−04	–

Table 23. The identified material properties of glass epoxy composite plate.

	Run 1	Run 2	Run 3	Run 4	Run 5	Average
E x (GPa)	44.3874	44.2024	44.0750	44.2620	44.5527	44.2959
E y (GPa)	17.0370	16.9951	17.1254	17.0207	17.1800	17.0716
G xy (GPa)	7.1892	7.2141	7.2229	7.2068	7.1833	7.2033
v xy	0.3100	0.3200	0.2856	0.3050	0.2901	0.3020
Obj. Fun.	3.3069–05	2.0609E−05	7.1041E−05	2.7600E−05	8.1624E−05	–

Table 24. The identified material properties of aluminium composite panel (ACP).

	Run 1	Run 2	Run 3	Run 4	Run 5	Average
E x (GPa)	12.5909	12.6524	12.6462	12.6805	12.3837	12.5907
E y (GPa)	10.3019	9.5653	10.0085	10.0935	9.7639	9.9466
G xy (GPa)	3.5114	3.6165	3.4915	3.4892	3.5649	3.5347
v xy (GPa)	0.3275	0.3126	0.3175	0.3352	0.3542	0.3294
Obj. Fun.	1.0647E−03	1.0813E−03	9.3888E−04	9.1873E−04	9.4841E−04	–

Table 25. Comparison of measured and evaluated natural frequencies of carbon epoxy composite plate, glass epoxy composite plate and aluminium composite panel (ACP) with F-F-F-F boundary condition.

Mode	Carbon epoxy			Glass epoxy			ACP
	Exp (Hz) [35]	Num (Hz)	Error (%)	Exp (Hz) [37]	Num (Hz)	Error (%)	Exp (Hz)	Num (Hz)	Error (%)
1	97.0	97.3589	0.3700	72.3	72.2547	0.0627	52.7	53.1642	0.8808
2	123.8	123.7941	0.0047	97.6	97.5105	0.0917	82.1	83.2382	1.3864
3	237.4	234.4539	1.2410	176.7	176.8566	0.0886	107.0	108.2136	1.1342
4	341.0	341.2816	0.0826	197.4	196.9956	0.2049	142.0	139.8169	1.5374
5	458.0	453.0001	1.0917	238.3	239.0003	0.2939	149.0	146.4892	1.6851
6	502.6	504.1175	0.3019	277.5	278.1660	0.2400	–	–	–
7	541.2	541.0719	0.0237	348.0	348.9600	0.2759	–	–	–
8	653.0	651.4398	0.2389	–	–	–	–	–	–

Notes: Exp: Experimental natural frequencies; Num: Numerical natural frequencies; Error: Absolute percentage error.

The credibility of the evaluated results should be validated with the results retrieved from an independent test. In the present study, destructive test results obtained from Riga Technical University (RTU), Israel Aircraft Industry (IAI) and German Aerospace Centre (DLR) are used as comparison references to validate the identified elastic properties of carbon epoxy composite plate, as shown in Table 26. It shows satisfactory agreement between the evaluated and reference results. For glass epoxy composite plate, non-destructive test results obtained from past literatures are used to validate the present evaluated results, as demonstrated in Table 27. The identified results can be seen coinciding considerably well with the reference results. In view of aluminium composite panel (ACP), the evaluated material properties show considerably good agreement with the benchmark material properties obtained from calculation using established formula, as presented in Table 28.

Table 26. Comparison of carbon epoxy composite plate elastic constants obtained from present test and destructive test references.

Elastic constant	Present	Destructive test references [35]
		RTU	IAI	DLR
E x (GPa)	170.5685	176	165	192
E y (GPa)	10.3249	8.9	9.2	10.6
G xy (GPa)	6.0626	5.2	5.4	6.1
v xy	0.3047	0.34	–	0.31

Notes: RTU: Riga Technical University; IAI: Israel Aircraft Industry; DLR: German Aerospace Centre.

Table 27. Comparison of glass epoxy composite plate elastic constants obtained from present test and non-destructive test references.

Elastic constant	Present	Non-destructive test references
		Rikards [36]	Arauio [38]	Sol. [39]
E x (GPa)	44.2024	44.2	43.6	45.6
E y (GPa)	16.9951	17.7	17.0	16.9
G xy (GPa)	7.2141	7.2	7.3	6.9
v xy	0.3200	0.195	0.350	0.220

Table 28. Comparison of aluminium composite panel (ACP) elastic constants obtained from present test and calculation.

	Present	Benchmark
E x (GPa)	12.6805	10.3097
E y (GPa)	10.0935	10.3097
G xy (GPa)	3.4892	3.8465
v xy	0.3352	0.3656

7. Conclusions

The paper presents a novel method, known as GA-ACO-PSO hybrid algorithm to effectively determine the elastic properties of various types of composite plates. The primary motivation of integrating GA, ACO and PSO is to minimize each other’s weaknesses and in the meantime, to promote respective strengths and therefore, transforming into a better hybrid algorithm. The feasibility of the proposed algorithm has been validated using published literatures as references. The use of three different types of boundary conditions during validation stage is investigated in order to study the reliability of the proposed algorithm. The proposed hybrid algorithm has been proven to be more reliable and effective in the aspects of accuracy, repeatability and convergence rate compared to GA, ACO and PSO. The robustness of the proposed algorithm to uncertainties (i.e. 1, 5 and 10% level of error) has also been successfully verified. The verified proposed algorithm has as well been used to identify the elastic properties of carbon epoxy composite plate, glass epoxy composite plate and aluminium composite panel (ACP). The outcomes turn out to be sufficiently good with absolute percentage errors between measured and evaluated natural frequencies not more than 2%. Therefore, the proposed hybrid algorithm is highly recommended for the identification of elastic properties of composite plates due to its better accuracy, repeatability, convergence rate as well as its robustness to uncertainties.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by Postgraduate Research Fund [grant number PG009-2015A]; Fundamental Research Grant Scheme [grant number FP010-2014A]; Advanced Shock and Vibration Research (ASVR) Group of University of Malaya; and other project collaborators.