# Memetic Computing Applied to the Design of Composite Materials and Structures

Authors: Jose Ignacio Pelaez; Jose Antonio Gomez-Ruiz; Jorge Veintimilla; Gustavo Vaccaro; Patricia Witt
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Presently, there exists an important need for lighter and more resistant structures, with reduced manufacturing costs. Laminated polymers are materials which respond to these new demands. Main difficulties of the design process of a composite laminate include the necessity to design both the geometry of the element and the material configuration itself and, therefore, the possibilities of creating composite materials are almost unlimited. Many techniques, ranging from linear programming or finite elements to computational intelligence, have been used to solve this type of problems. The aim of this work is to show that more effective and dynamic methods to solve this type of problems are obtained by using certain techniques based on systematic exploitation of knowledge of the problem, together with the combination of metaheuristics based on population as well as on local search. With this objective, a memetic algorithm has been designed and compared with the main heuristics used in the design of laminated polymers in different scenarios. All solutions obtained have been validated by the ANSYS® software package.


---

## Body

## 1. Introduction

A composite material is formed by the aggregation of two or more distinct materials to form a new one with enhanced properties: an agglomerate material known as the matrix and reinforcement materials that may be made up of continuous fibers, short fibers, or particles [1, 2]. Well-designed material adopts the best properties of its constituents and even some that none of these possess. The aim of composite materials design is to generate new cheaper materials with improved strength and lightness. Not all of them can be improved simultaneously; therefore, the design objective is to obtain a new material that offers the best possible adaptation to the required specifications.The exceptional strength and lightness of these materials have led to the development of a vast number of applications, particularly in the aeronautical and space industries due to the economic significance of these properties. Figure1 shows the markets for composite materials.Figure 1
Markets for composite materials.Laminate is a particularly important type of composite material made up of laminas of the same composition with unidirectional reinforcement fibers, stacked, and bound together by the same material that forms the matrix, but with distinct fiber orientation. Within this category, symmetric laminates are worthy of special attention, with both geometrical and structural symmetry relative to the mid-plane.One of the main difficulties of the design process of a laminate is to design both the geometry of the element and the material configuration itself as to best exploit the qualities of the constituent materials. Designing process must also evaluate the deterioration of the laminate properties over time due to stress, which can lead to unanticipated behaviour (cracks) or failure of the structural element in question.Synthesis and analysis have been traditionally carried out using empirical knowledge based methods [3]. This is partly because the number of possible combinations of composites is almost unlimited and also because characterization by experimentation is very expensive.Since the 1990s, different design systems which aim to overcome these limitations have been proposed. These proposals have involved approaches ranging from traditional techniques such as classical nonlinear optimization procedures combined with finite element modelling [4, 5], through generic task methods and case-based reasoning [6], to modern artificial intelligence techniques [7–16].Memetic Computing or Memetic Algorithms (MA) have proven to be efficient at numerous situations [17–21]. Figure 2 shows the scheme for a typical MA [22]. Local search is performed between each generation, in addition to the techniques used by Genetic Algorithms to explore the search space, namely, recombination/crossover and mutation. It is performed to improve the fitness of the population (in a localized region of the solution space) so that the next generation has “better” genes from its parents.Figure 2
Scheme of the memetic algorithm.fth denotes a fixed threshold value used for the stopping criterion.The aim of this work is to show that more effective and dynamic methods to solve problems of laminated polymers design are obtained by using certain techniques based on systematic exploitation of knowledge of the problem, together with the combination of metaheuristics based on population as well as on local search. With this objective, a MA has been designed and thoroughly compared with the main heuristics used in the design of laminated polymers in different scenarios. This work has been organized as follows: in Section2, a studio of the main metaheuristics which have been used in the design of composite laminates is performed; in Section 3, a model for the design and optimization of symmetric laminates based on a MA is presented; in Section 4, this model is subjected to an exhaustive analysis and is compared with other methods to prove its usefulness; in Section 5, the MA-based model is applied to specific cases. Finally, conclusions are presented in Section 6.


## 2. State of the Art of Composite Materials Design

### 2.1. Genetic Algorithms

Genetic Algorithms (GA) have been the most popular method in the design of laminated polymers and optimizing piling sequence [23]. Callahan and Weeks [24], Nagendra et al. [25], Le Riche and Haftka [26], and Ball et al. [27] have been the first to adopt and use GA for the design of piling sequences of composite laminated materials. GA has also been used in problems with different objective functions, such as strength [28, 29], buckling loads [9, 28, 30–34], dimensional stability [35], strain energy absorption [36], weight (either as a restriction or as an objective to minimize) [37, 38], bending/torsion connection, stiffness [36, 39], basic frequencies [34, 40–42], distortion [43], or finding laminate reference parameters [44].GA have also been applied in the design of a variety of composite structures ranging from simple rectangular plates to complex geometrical sheets, such as sandwich panels [45], rigid sheets [46], bolt joints [47], and laminated cylindrical panels [34]. Similarly, GA have been combined with finite elements packages which analyse the response to tension and deformation of the composite material structure [43, 48, 49]. A combination of methods is sometimes used; for example, Park et al. [50] uses an approximate memory combined with a permutation operator and with the local learning/random mixture, in order to reduce function number and improve rate of convergence.Some of the main problems regarding GA are deep computational necessity and premature convergence. In order to solve these drawbacks, increase the rate of convergence, reduce risk of premature convergence, and lessen function evaluation time, several modifications have been proposed, including the use of parallel computation [51, 52]; optimization of several levels (thick and thin level codification) [51]; introduction of problem-dependable operators [26]; layer adding or layer deleting, permutation, and interlaminar change [9, 46]; generalized elitist and mutation of thickness of laminate/material/angle of fiber [49]; recovery of previously evaluated solutions [10, 50]; using approximation methods for the evaluation function [33, 53, 54]; training artificial neural network [34, 42]; initial blood-related population; or an aging hierarchic structure [55].Sargent et al. [56] compared GA with some raging algorithms (i.e., random search, raging search, and simulated annealing) and observed that GA obtained better solutions than raging search, but in some cases they were unable to find a solution. Sivakumar et al. [57] compared David, Fletcher, and Powell (DFP)’s Quasi-Newton Method and GA, applied to reduce weight of a laminated composite limited by its basic frequencies. It was concluded that DFP converged in a reduced number of iterations when restriction number was small. However, finding a feasible number was difficult when restrictions number increased. Also considering that DFP could not handle discreet variables, it was concluded that “GA was a better tool to optimize laminated composites.”Even though GA has been widely used in optimization of piling sequence, an important flaw is its low rate of convergence. GA is an evolutionary algorithm based on population and might require several generations before converging into a solution [58]. Each generation consists of a great number of function evaluations; therefore, it could require plenty of computational time besides being very expensive.


### 2.2. Simulated Annealing

After Genetic Algorithms, simulated annealing (SA) is the second most popular method for the optimization of piling sequence in laminated composite materials [59–62].The main problem of this technique is the generation of a sequence of points that converges into a nonoptimal solution. In order to solve this lack, some modifications have been proposed, such as an increase in the probability of sample points far from the present point [63] or the use of a set of points at a time instead of only one [64]. In order to increase the rate of convergence, Genovese et al. [65] proposed a two-level algorithm, including a “global annealing” in which all the design variables were disrupted simultaneously and a “local annealing” in which only one design variable was disrupted at a time. The local annealing was performed after each iteration of the global annealing in an attempt to improve the testing point locally. It was found that its speed of convergence was greater than the one-level simulated annealing and comparable with the optimization method based on gradient implemented by the Sequential Linear Programming (SLP) method.SA is a good choice for the general case of optimal selection of laminate; however, it cannot be programmed to take benefit of the advantages of specific properties of a certain problem. GA is, in this respect, more flexible, but it frequently takes more computational time [56, 61, 66]. It is not easy to generalize the previous conclusion, since there are other researches, such as Rama Mohan Rao and Shyju [58], which shows that SA had better computational efficiency and was better at finding a solution to other combinatorial problems.


### 2.3. Tabu Search

Tabu Search (TS) was implemented by Pai et al. [14] for the optimization of piling sequence of laminates subjected to bulging and resistance requests as well as for matrix rupture. Results were compared with GA and they showed a comparative solution, but computational time depended on the case.It has been used combined with other techniques, such as the SA. Rama Mohan Rao and Arvind [67] added a TS in the SA obtaining a method called tabu embedded in simulated annealing (TSA). Optimization of piling sequence of laminated composites was solved by the TSA, whilst restrictions are administrated through a correction strategy. TSA was faster than classic SA although more memory is required.


### 2.4. Variable Neighbourhood Search

Variable neighbourhood search (VNS) has been recently applied to laminated polymers design. Corz et al. [68] propose and compare an algorithm based on a variable environments search that allows the design of the geometry and composition of symmetric laminated composite materials. In this work, the implementation of the algorithm is shown in detail: objective function, encoding, fitness, neighbourhood structures, and local search, as well as the data used in the examples: volume fraction, efforts, and coefficients. The proposed model is compared with other techniques such as GA, SA, and TS, showing more efficient results.


## 2.1. Genetic Algorithms

Genetic Algorithms (GA) have been the most popular method in the design of laminated polymers and optimizing piling sequence [23]. Callahan and Weeks [24], Nagendra et al. [25], Le Riche and Haftka [26], and Ball et al. [27] have been the first to adopt and use GA for the design of piling sequences of composite laminated materials. GA has also been used in problems with different objective functions, such as strength [28, 29], buckling loads [9, 28, 30–34], dimensional stability [35], strain energy absorption [36], weight (either as a restriction or as an objective to minimize) [37, 38], bending/torsion connection, stiffness [36, 39], basic frequencies [34, 40–42], distortion [43], or finding laminate reference parameters [44].GA have also been applied in the design of a variety of composite structures ranging from simple rectangular plates to complex geometrical sheets, such as sandwich panels [45], rigid sheets [46], bolt joints [47], and laminated cylindrical panels [34]. Similarly, GA have been combined with finite elements packages which analyse the response to tension and deformation of the composite material structure [43, 48, 49]. A combination of methods is sometimes used; for example, Park et al. [50] uses an approximate memory combined with a permutation operator and with the local learning/random mixture, in order to reduce function number and improve rate of convergence.Some of the main problems regarding GA are deep computational necessity and premature convergence. In order to solve these drawbacks, increase the rate of convergence, reduce risk of premature convergence, and lessen function evaluation time, several modifications have been proposed, including the use of parallel computation [51, 52]; optimization of several levels (thick and thin level codification) [51]; introduction of problem-dependable operators [26]; layer adding or layer deleting, permutation, and interlaminar change [9, 46]; generalized elitist and mutation of thickness of laminate/material/angle of fiber [49]; recovery of previously evaluated solutions [10, 50]; using approximation methods for the evaluation function [33, 53, 54]; training artificial neural network [34, 42]; initial blood-related population; or an aging hierarchic structure [55].Sargent et al. [56] compared GA with some raging algorithms (i.e., random search, raging search, and simulated annealing) and observed that GA obtained better solutions than raging search, but in some cases they were unable to find a solution. Sivakumar et al. [57] compared David, Fletcher, and Powell (DFP)’s Quasi-Newton Method and GA, applied to reduce weight of a laminated composite limited by its basic frequencies. It was concluded that DFP converged in a reduced number of iterations when restriction number was small. However, finding a feasible number was difficult when restrictions number increased. Also considering that DFP could not handle discreet variables, it was concluded that “GA was a better tool to optimize laminated composites.”Even though GA has been widely used in optimization of piling sequence, an important flaw is its low rate of convergence. GA is an evolutionary algorithm based on population and might require several generations before converging into a solution [58]. Each generation consists of a great number of function evaluations; therefore, it could require plenty of computational time besides being very expensive.


## 2.2. Simulated Annealing

After Genetic Algorithms, simulated annealing (SA) is the second most popular method for the optimization of piling sequence in laminated composite materials [59–62].The main problem of this technique is the generation of a sequence of points that converges into a nonoptimal solution. In order to solve this lack, some modifications have been proposed, such as an increase in the probability of sample points far from the present point [63] or the use of a set of points at a time instead of only one [64]. In order to increase the rate of convergence, Genovese et al. [65] proposed a two-level algorithm, including a “global annealing” in which all the design variables were disrupted simultaneously and a “local annealing” in which only one design variable was disrupted at a time. The local annealing was performed after each iteration of the global annealing in an attempt to improve the testing point locally. It was found that its speed of convergence was greater than the one-level simulated annealing and comparable with the optimization method based on gradient implemented by the Sequential Linear Programming (SLP) method.SA is a good choice for the general case of optimal selection of laminate; however, it cannot be programmed to take benefit of the advantages of specific properties of a certain problem. GA is, in this respect, more flexible, but it frequently takes more computational time [56, 61, 66]. It is not easy to generalize the previous conclusion, since there are other researches, such as Rama Mohan Rao and Shyju [58], which shows that SA had better computational efficiency and was better at finding a solution to other combinatorial problems.


## 2.3. Tabu Search

Tabu Search (TS) was implemented by Pai et al. [14] for the optimization of piling sequence of laminates subjected to bulging and resistance requests as well as for matrix rupture. Results were compared with GA and they showed a comparative solution, but computational time depended on the case.It has been used combined with other techniques, such as the SA. Rama Mohan Rao and Arvind [67] added a TS in the SA obtaining a method called tabu embedded in simulated annealing (TSA). Optimization of piling sequence of laminated composites was solved by the TSA, whilst restrictions are administrated through a correction strategy. TSA was faster than classic SA although more memory is required.


## 2.4. Variable Neighbourhood Search

Variable neighbourhood search (VNS) has been recently applied to laminated polymers design. Corz et al. [68] propose and compare an algorithm based on a variable environments search that allows the design of the geometry and composition of symmetric laminated composite materials. In this work, the implementation of the algorithm is shown in detail: objective function, encoding, fitness, neighbourhood structures, and local search, as well as the data used in the examples: volume fraction, efforts, and coefficients. The proposed model is compared with other techniques such as GA, SA, and TS, showing more efficient results.


## 3. A Memetic Algorithm-Based Model for the Designing of Composite Materials

In this section, the optimization problem, an encoding scheme for the representation of the problem, a fitness function to evaluate the feasibility of solutions, and the reproductive and local search operators used are presented.


### 3.1. Optimization Problem

The optimization problem can be formulated as follows: find the material (fiberF and matrix M), the volume fraction, VF, and the laminate stacking sequence, θ1,…,θk,…,θn, with the purpose of maximizing the utilization of the laminate, thus achieving the lowest number of laminas n. The set of design variables is expressed as the vector s~=(F,M,VF,θ1,…,θk,…,θn). The degree of utilization of the laminate material is defined as a positive real number that is obtained from the fitness function (FF), defined in the Section 3.2, to the s~ vector. Therefore, the optimization problem and its corresponding constrains can be defined as follows:(1)maxFFs~,where n≥2 (number of laminas);VF∈0.3,0.7 in 0.1 intervals (volume fraction);θi∈-80°,90° in 10° intervals (orientation of the laminas).


### 3.2. Encoding

The structure of the laminate will be represented by an agent formed by four information units (IU), which encode the fiber, the matrix, the volume fraction, and the laminate geometry (see Figure3). Fitness, that shows how well adapted the solution is to its environment, is added to the structure.Figure 3
Structure of the laminate.Each of the four IU that define the laminate represents a different characteristic of its composition: fiber and matrix are represented by an integer number; volume fraction is represented by a real number; and geometry is represented by a sequence of integers that represent the direction of fibers. The values used for the different IU are shown in Table1.Table 1
Encoding used for a laminate.
IU
Values
Fiber
A set of eleven fibres are used, each represented by an integer number
1
2
3
4
5
6
7
8
9
10
11
E-Glass
S-Glass
AS-1
AS-4
IM-7
P-100
T-40
T-300
Boron
Kevlar-49
Kevlar-149
Matrix
A set of five matrices are used, each represented by an integer number
1
2
3
4
5
Peek
PPS
Polyester
Epoxy
Polyamide
Volume fraction
A real number from the set0.3,0.4,0.5,0.6,0.7
Geometry of the laminate
A sequence of integers, one for each lamina, in the range-80°,+90°, in 10° intervalsTo represent the geometry of the laminate, we use the following notation: [50/−10/0/60/−20]s. It is a symmetric laminate whose outer (first) lamina has an orientation of 50°, the second one (next) of −10°, and so on. A simplified representation of a symmetric laminate is shown in Figure 4.Figure 4
Symmetric laminate.Different failure criteria can be used to determine whether a lamina can withstand specific stress conditions without breaking [1]. In this article, the Tsai-Wu failure criteria [69] have been selected. In order to determine if a lamina breaks, for each lamina i of a laminate, its Tsai-Wu coefficient Pki is defined as(2)Pki=1X-1X′·σ1+1Y-1Y′·σ2+1X·X′·σ12+1Y·Y′·σ22-1X·X′·Y·Y′·σ1·σ2+1S2·τ122,where X,X′,Y,Y′, and S are, respectively, the ultimate tensile and compressive strength in the fiber direction, the ultimate tensile and compressive transversal strength, and the shear strength. σ1 and σ2 are the stress resistance coefficients for the x- and y-axes, respectively, and τ12 is the stress coefficient of an orthotropic lamina under plane stress conditions. Applying the Tsai-Wu failure criteria [69], the lamina does not break if Pki≤1, ∀i, and breaks otherwise. Furthermore, if Pki≅1, the lamina i is working at full capacity. If a lamina is broken, then the entire laminate is discarded.


### 3.3. Fitness

The fitness function considers economic and safety criteria, so the alignment of stresses along the direction of the fibers is improved and the following are penalized: volume fraction and high laminate thickness; stacking more than four consecutive laminas with the same fiber orientation; the distribution of stress in directions other than that of the fibers; and lamina breakage. A higher value indicates a better solution. The fitness function,FF, is defined as(3)FF=10αP1βCLA+1·VF·n·eγ·R+1δ·P2+1·P12+1,where the coefficients are described in Table 2. The multiplicative factor 10∝ is to maintain fitness values between manageable values. α has been empirically determined as −28.Table 2
Fitness equation coefficients.
Coefficient
Description
P
1
Longitudinal coefficient along the direction of the fibres. It is a real value indicating the relationship between longitudinal stress in the direction of the fibres and the breakage stress in that direction. As the direction of the fibres is that of most resistance, it is of interest that this ratio is as high as possible. For this reason, it is in the numerator.β has been empirically determined as 2.P1=∑iσ1i/Xi∨σ1i/X′i where the stress σ1i is either tensile or compressive and the sum applies to the laminas that do not break.
P
2
Longitudinal coefficient perpendicular to the direction of the fibres. It is a real value that indicates the relationship between the longitudinal stress perpendicular to the fibres and the breakage stress in that direction. As the direction which is perpendicular to the fibres is that of least resistance, it is of interest that this ratio is as small as possible. For this reason it is in the denominator.P2=∑iσ2i/Yi∨σ2i/Y′i where the stress σ2i is either tensile or compressive and the sum applies to all the laminas.
P
12
Shear coefficient. It is a real value that indicates the relationship between the shear stress in the plane and the maximum stress that can be tolerated in that plane. As the lamina has low resistance to these stresses, it is of interest that this ratio is as small as possible. For this reason it is in the denominator. The sum applies to all the laminas.P12=∑iτ12i/Si
C
L
A
Coefficient that indicates whether the number of adjacent layers with the same orientation is less than or equal to fourCLA=0 or greater than four CLA=1010 It is in the denominator so as to heavily penalize this second possibility. Its values have been empirically determined.
V
F
Fiber volume fraction of the laminate. As this value should tend towards small values that are more advantageous economically, it is placed in the denominator where small values will have a positive effect.
n
·
e
Product of the number of laminas and its thickness (expressed in metres).γ has been empirically determined as 12.
R
It indicates the number of layers that breakPki>1, and is in the denominator in order to penalize such event. δ has been empirically determined as 4.


### 3.4. Reproductive Operator

Crossover operator is different for each IU In case of fiber, matrix, and volume fraction IU, a classic one-point crossover operator adapted to the different possible cases is used. However, for geometry IU, the classic crossover operator is not used because it is very static in changing the number of layers. To avoid this problem, a dynamic crossover operator with a different crossover point for each parent is proposed. Such crossover points are determined by generating two integer numbers from a continuous uniform distribution. The first number lies between 0 and the number of laminas of the first laminate,l1, and the second between 1 and the number of laminas of the second laminate, l2, plus 1. Then, the parent laminates are split into two parts, with the child joining together the left part of the first parent and the right part of the second parent being formed. Depending on the location of these crossover points, p1 and p2 can be obtained different cases shown in Table 3.Table 3
Different cases for crossover operator.
Case
Scheme
(1) This is the general case. The child is formed by joining the1⋯p1 laminas of the first parent, with the p2⋯l2 laminas of the second one.If the length of the child exceeds the maximum allowed, it is cut to achieve this ultimate size.
(2) If the first crossover point,p1, is zero, and the second crossover point, p2, is not l2+1, the child is formed with the p2⋯l2 laminas of the second parent.
(3) If the first crossover point,p1, is not zero, and the second crossover point, p2, is l2+1, the child is formed with the 1⋯p1 laminas of the first parent.
(4) If the first crossover point,p1, is zero, and the second crossover point, p2, is l2+1, the child is formed with the first lamina of the first parent.


### 3.5. Local Search Operators

A set of local search operators are defined, following previous applications of a variable neighbourhood search-based model [68]. These operators allow a lamina to be added or removed, to change its orientation, type of matrix, and type of fiber or volume fraction. Two new operators have been defined in order to better exploit the region being explored in a given time. The operators proposed in the model are summarized in Table 4.Table 4
Local search operators.
Operator
Description
Function
Number 1
It operates on the geometry IU of the laminate, varying the orientation of the laminas, moving each of them 10° closer to 0°. If the orientation of a lamina is 0°, it is left unchanged.
To obtain solutions with improved longitudinal stress tolerance along thex-axis better.
Number 2
It operates on the geometry IU of the laminate, varying the orientation of the laminas, moving each of them 10° closer to 90°. If the orientation of a lamina is 90°, it is left unchanged, and if it is −80° it is changed to 90°.
To obtain solutions with improved longitudinal stress tolerance along they-axis better.
Number 3
It operates on the geometry IU of the laminate, varying the orientation of the laminas, moving each of them 10° closer to +45° or −45°, depending on whether it is positive or negative. If the orientation of a lamina is 40°, 50°, −40° or −50°, it is left unchanged, and if it is 0° it is changed to 10°.
To obtain solutions with improved shear stress tolerance.
Number 4
It operates on the geometry IU of the laminate, adding an individual lamina with random orientation and positioning. If the number of laminas is already the maximum possible, it is left unchanged.
To explore into the direction of the greatest number of laminas
Number 5
It operates on the geometry I.U. of the laminate, removing an individual lamina, chosen at random. If there is only one lamina, it is left unchanged.
To explore into the direction of the least number of laminas
Number 6
It operates on the geometry IU of the laminate, changing the orientation of a lamina, chosen at random and relocated in a new position, also chosen at random.
To explore into the direction of new orientations.
Number 7
It operates on the fiber IU, changing it by a distinct one, chosen at random, maintaining the matrix, the volume fraction, the number of laminas and their orientation constant.
To explore into the subspace solution of the fibres alone.
Number 8
It operates on the matrix I.U., changing it by a distinct one, chosen at random, maintaining the fiber, the volume fraction, the number of laminas and their orientation constant.
To explore into the subspace solution of the matrices alone.
Number 9
It operates on the volume fraction IU, changing its value for a distinct one, chosen at random, within the acceptable limits, maintaining the fiber, the matrix, the number of laminas and their orientation constant.
To explore into the subspace solution of the volume fraction alone.
Number 10
It operates on the geometry IU of the laminate, adding an individual lamina with orientation of −10° from the last positioned lamina. If the number of laminas is already the maximum possible, it is left unchanged.
To explore into the greatest number of laminas.
Number 11
It operates on the geometry IU of the laminate, removing an individual lamina, the closest to the last removed lamina, with different angular orientation. If there is only one lamina, it is left unchanged.
To explore into the smallest number of laminas.


### 3.6. Implementation

The MA proposed model is compared with other heuristics in order to verify its advantages. These heuristics (detailed in Section2) are GA, VNS, SA, and TS. All heuristics use the same representation of the problem proposed in this paper. Experiments were performed in the Picasso Computer at University of Malaga (512 IBM PowerPC 970 processors, 2994.04 GFlops).In Tables5 and 6, characteristics of the materials (fibers and matrices) used in the tests for the laminates design are shown.Table 5
Characteristics of fibers used in the design process.
Material
Module elasticE[GPa]
Resistance tractionσf[MPa]
Coefficient Poissonυ
Glass fiber
E-Glass
72.40
3450.00
0.20
S_Glass
86.90
4300.00
0.22
Carbon fiber
AS-1
228.00
3100.00
0.20
AS-4
248.00
4070.00
0.20
IM-7
301.00
5310.00
0.20
P-100
758.00
2410.00
0.20
T-40
290.00
5650.00
0.20
T-300
231.00
3650.00
0.20
Boron fiber
Boron
393.00
3100.00
0.20
Aramid fiber
Kevlar 49
131.00
3620.00
0.35
Kevlar 149
179.00
3450.00
0.35Table 6
Characteristics of matrices used in the design process.
Material
Module elasticE[GPa]
Resistance tractionσf[MPa]
Coefficient Poissonυ
Thermoplastics
Peek
3.24
100.00
0.40
PPS
3.30
82.70
0.33
Thermostables
Polyester
3.00
50.00
0.30
Epoxy resin
4.60
58.60
0.36
Polyamide
3.50
103.00
0.35


## 3.1. Optimization Problem

The optimization problem can be formulated as follows: find the material (fiberF and matrix M), the volume fraction, VF, and the laminate stacking sequence, θ1,…,θk,…,θn, with the purpose of maximizing the utilization of the laminate, thus achieving the lowest number of laminas n. The set of design variables is expressed as the vector s~=(F,M,VF,θ1,…,θk,…,θn). The degree of utilization of the laminate material is defined as a positive real number that is obtained from the fitness function (FF), defined in the Section 3.2, to the s~ vector. Therefore, the optimization problem and its corresponding constrains can be defined as follows:(1)maxFFs~,where n≥2 (number of laminas);VF∈0.3,0.7 in 0.1 intervals (volume fraction);θi∈-80°,90° in 10° intervals (orientation of the laminas).


## 3.2. Encoding

The structure of the laminate will be represented by an agent formed by four information units (IU), which encode the fiber, the matrix, the volume fraction, and the laminate geometry (see Figure3). Fitness, that shows how well adapted the solution is to its environment, is added to the structure.Figure 3
Structure of the laminate.Each of the four IU that define the laminate represents a different characteristic of its composition: fiber and matrix are represented by an integer number; volume fraction is represented by a real number; and geometry is represented by a sequence of integers that represent the direction of fibers. The values used for the different IU are shown in Table1.Table 1
Encoding used for a laminate.
IU
Values
Fiber
A set of eleven fibres are used, each represented by an integer number
1
2
3
4
5
6
7
8
9
10
11
E-Glass
S-Glass
AS-1
AS-4
IM-7
P-100
T-40
T-300
Boron
Kevlar-49
Kevlar-149
Matrix
A set of five matrices are used, each represented by an integer number
1
2
3
4
5
Peek
PPS
Polyester
Epoxy
Polyamide
Volume fraction
A real number from the set0.3,0.4,0.5,0.6,0.7
Geometry of the laminate
A sequence of integers, one for each lamina, in the range-80°,+90°, in 10° intervalsTo represent the geometry of the laminate, we use the following notation: [50/−10/0/60/−20]s. It is a symmetric laminate whose outer (first) lamina has an orientation of 50°, the second one (next) of −10°, and so on. A simplified representation of a symmetric laminate is shown in Figure 4.Figure 4
Symmetric laminate.Different failure criteria can be used to determine whether a lamina can withstand specific stress conditions without breaking [1]. In this article, the Tsai-Wu failure criteria [69] have been selected. In order to determine if a lamina breaks, for each lamina i of a laminate, its Tsai-Wu coefficient Pki is defined as(2)Pki=1X-1X′·σ1+1Y-1Y′·σ2+1X·X′·σ12+1Y·Y′·σ22-1X·X′·Y·Y′·σ1·σ2+1S2·τ122,where X,X′,Y,Y′, and S are, respectively, the ultimate tensile and compressive strength in the fiber direction, the ultimate tensile and compressive transversal strength, and the shear strength. σ1 and σ2 are the stress resistance coefficients for the x- and y-axes, respectively, and τ12 is the stress coefficient of an orthotropic lamina under plane stress conditions. Applying the Tsai-Wu failure criteria [69], the lamina does not break if Pki≤1, ∀i, and breaks otherwise. Furthermore, if Pki≅1, the lamina i is working at full capacity. If a lamina is broken, then the entire laminate is discarded.


## 3.3. Fitness

The fitness function considers economic and safety criteria, so the alignment of stresses along the direction of the fibers is improved and the following are penalized: volume fraction and high laminate thickness; stacking more than four consecutive laminas with the same fiber orientation; the distribution of stress in directions other than that of the fibers; and lamina breakage. A higher value indicates a better solution. The fitness function,FF, is defined as(3)FF=10αP1βCLA+1·VF·n·eγ·R+1δ·P2+1·P12+1,where the coefficients are described in Table 2. The multiplicative factor 10∝ is to maintain fitness values between manageable values. α has been empirically determined as −28.Table 2
Fitness equation coefficients.
Coefficient
Description
P
1
Longitudinal coefficient along the direction of the fibres. It is a real value indicating the relationship between longitudinal stress in the direction of the fibres and the breakage stress in that direction. As the direction of the fibres is that of most resistance, it is of interest that this ratio is as high as possible. For this reason, it is in the numerator.β has been empirically determined as 2.P1=∑iσ1i/Xi∨σ1i/X′i where the stress σ1i is either tensile or compressive and the sum applies to the laminas that do not break.
P
2
Longitudinal coefficient perpendicular to the direction of the fibres. It is a real value that indicates the relationship between the longitudinal stress perpendicular to the fibres and the breakage stress in that direction. As the direction which is perpendicular to the fibres is that of least resistance, it is of interest that this ratio is as small as possible. For this reason it is in the denominator.P2=∑iσ2i/Yi∨σ2i/Y′i where the stress σ2i is either tensile or compressive and the sum applies to all the laminas.
P
12
Shear coefficient. It is a real value that indicates the relationship between the shear stress in the plane and the maximum stress that can be tolerated in that plane. As the lamina has low resistance to these stresses, it is of interest that this ratio is as small as possible. For this reason it is in the denominator. The sum applies to all the laminas.P12=∑iτ12i/Si
C
L
A
Coefficient that indicates whether the number of adjacent layers with the same orientation is less than or equal to fourCLA=0 or greater than four CLA=1010 It is in the denominator so as to heavily penalize this second possibility. Its values have been empirically determined.
V
F
Fiber volume fraction of the laminate. As this value should tend towards small values that are more advantageous economically, it is placed in the denominator where small values will have a positive effect.
n
·
e
Product of the number of laminas and its thickness (expressed in metres).γ has been empirically determined as 12.
R
It indicates the number of layers that breakPki>1, and is in the denominator in order to penalize such event. δ has been empirically determined as 4.


## 3.4. Reproductive Operator

Crossover operator is different for each IU In case of fiber, matrix, and volume fraction IU, a classic one-point crossover operator adapted to the different possible cases is used. However, for geometry IU, the classic crossover operator is not used because it is very static in changing the number of layers. To avoid this problem, a dynamic crossover operator with a different crossover point for each parent is proposed. Such crossover points are determined by generating two integer numbers from a continuous uniform distribution. The first number lies between 0 and the number of laminas of the first laminate,l1, and the second between 1 and the number of laminas of the second laminate, l2, plus 1. Then, the parent laminates are split into two parts, with the child joining together the left part of the first parent and the right part of the second parent being formed. Depending on the location of these crossover points, p1 and p2 can be obtained different cases shown in Table 3.Table 3
Different cases for crossover operator.
Case
Scheme
(1) This is the general case. The child is formed by joining the1⋯p1 laminas of the first parent, with the p2⋯l2 laminas of the second one.If the length of the child exceeds the maximum allowed, it is cut to achieve this ultimate size.
(2) If the first crossover point,p1, is zero, and the second crossover point, p2, is not l2+1, the child is formed with the p2⋯l2 laminas of the second parent.
(3) If the first crossover point,p1, is not zero, and the second crossover point, p2, is l2+1, the child is formed with the 1⋯p1 laminas of the first parent.
(4) If the first crossover point,p1, is zero, and the second crossover point, p2, is l2+1, the child is formed with the first lamina of the first parent.


## 3.5. Local Search Operators

A set of local search operators are defined, following previous applications of a variable neighbourhood search-based model [68]. These operators allow a lamina to be added or removed, to change its orientation, type of matrix, and type of fiber or volume fraction. Two new operators have been defined in order to better exploit the region being explored in a given time. The operators proposed in the model are summarized in Table 4.Table 4
Local search operators.
Operator
Description
Function
Number 1
It operates on the geometry IU of the laminate, varying the orientation of the laminas, moving each of them 10° closer to 0°. If the orientation of a lamina is 0°, it is left unchanged.
To obtain solutions with improved longitudinal stress tolerance along thex-axis better.
Number 2
It operates on the geometry IU of the laminate, varying the orientation of the laminas, moving each of them 10° closer to 90°. If the orientation of a lamina is 90°, it is left unchanged, and if it is −80° it is changed to 90°.
To obtain solutions with improved longitudinal stress tolerance along they-axis better.
Number 3
It operates on the geometry IU of the laminate, varying the orientation of the laminas, moving each of them 10° closer to +45° or −45°, depending on whether it is positive or negative. If the orientation of a lamina is 40°, 50°, −40° or −50°, it is left unchanged, and if it is 0° it is changed to 10°.
To obtain solutions with improved shear stress tolerance.
Number 4
It operates on the geometry IU of the laminate, adding an individual lamina with random orientation and positioning. If the number of laminas is already the maximum possible, it is left unchanged.
To explore into the direction of the greatest number of laminas
Number 5
It operates on the geometry I.U. of the laminate, removing an individual lamina, chosen at random. If there is only one lamina, it is left unchanged.
To explore into the direction of the least number of laminas
Number 6
It operates on the geometry IU of the laminate, changing the orientation of a lamina, chosen at random and relocated in a new position, also chosen at random.
To explore into the direction of new orientations.
Number 7
It operates on the fiber IU, changing it by a distinct one, chosen at random, maintaining the matrix, the volume fraction, the number of laminas and their orientation constant.
To explore into the subspace solution of the fibres alone.
Number 8
It operates on the matrix I.U., changing it by a distinct one, chosen at random, maintaining the fiber, the volume fraction, the number of laminas and their orientation constant.
To explore into the subspace solution of the matrices alone.
Number 9
It operates on the volume fraction IU, changing its value for a distinct one, chosen at random, within the acceptable limits, maintaining the fiber, the matrix, the number of laminas and their orientation constant.
To explore into the subspace solution of the volume fraction alone.
Number 10
It operates on the geometry IU of the laminate, adding an individual lamina with orientation of −10° from the last positioned lamina. If the number of laminas is already the maximum possible, it is left unchanged.
To explore into the greatest number of laminas.
Number 11
It operates on the geometry IU of the laminate, removing an individual lamina, the closest to the last removed lamina, with different angular orientation. If there is only one lamina, it is left unchanged.
To explore into the smallest number of laminas.


## 3.6. Implementation

The MA proposed model is compared with other heuristics in order to verify its advantages. These heuristics (detailed in Section2) are GA, VNS, SA, and TS. All heuristics use the same representation of the problem proposed in this paper. Experiments were performed in the Picasso Computer at University of Malaga (512 IBM PowerPC 970 processors, 2994.04 GFlops).In Tables5 and 6, characteristics of the materials (fibers and matrices) used in the tests for the laminates design are shown.Table 5
Characteristics of fibers used in the design process.
Material
Module elasticE[GPa]
Resistance tractionσf[MPa]
Coefficient Poissonυ
Glass fiber
E-Glass
72.40
3450.00
0.20
S_Glass
86.90
4300.00
0.22
Carbon fiber
AS-1
228.00
3100.00
0.20
AS-4
248.00
4070.00
0.20
IM-7
301.00
5310.00
0.20
P-100
758.00
2410.00
0.20
T-40
290.00
5650.00
0.20
T-300
231.00
3650.00
0.20
Boron fiber
Boron
393.00
3100.00
0.20
Aramid fiber
Kevlar 49
131.00
3620.00
0.35
Kevlar 149
179.00
3450.00
0.35Table 6
Characteristics of matrices used in the design process.
Material
Module elasticE[GPa]
Resistance tractionσf[MPa]
Coefficient Poissonυ
Thermoplastics
Peek
3.24
100.00
0.40
PPS
3.30
82.70
0.33
Thermostables
Polyester
3.00
50.00
0.30
Epoxy resin
4.60
58.60
0.36
Polyamide
3.50
103.00
0.35


## 4. Experiments and Results

The problem is to determine the composition and thickness of a laminated plate under a distributedNx and Ny loading. Table 7 presents the statistics for 21 different loading cases (n=21) and the 5 heuristics (k=5) indicated in Section 3.6, with 500 simulations each. For each model, the results shown correspond to the smallest number of laminas of the material obtained, NLmin, the number of times this smallest number is obtained, NL, and its average NL¯. Also, the best (maximum) fitness obtained, FFmax, and its average FF¯ are presented.Table 7
Statistics corresponding to the series of simulations. Biaxial loadingNx:Ny.
Loads
Algorithm
N
L
m
i
n
NL
N
L
¯
F
F
m
a
x
F
F
¯
N
x: 250.000N/m Ny: 250.000N/m
VNS
6
500
6
2.8594·108
2.7272·108
MA
6
500
6
2.8594·108
2.8594·108
AG
6
89
8.28
2.8594·108
6.969·107
SA
8
76
12.456
8.978·106
1.1579·105
TS
6
32
8.7120
2.7920·108
7.4150·106
N
x: 250.000N/m Ny: 500.000N/m
VNS
8
500
8
1.4029·107
1.2698·107
MA
8
500
8
1.4029·107
1.3822·107
AG
8
35
11.8
1.3145·107
1.2643·106
SA
8
10
16.544
5.6258·106
2.0485·104
TS
8
13
10.5680
1.0056·107
1.5439·105
N
x: 250.000N/m Ny: 750.000N/m
VNS
8
498
8.008
6.6161·106
6.3187·106
MA
8
500
8
6.6161·106
6.6044·106
AG
8
3
14.4733
6.5042·106
1.9654·105
SA
10
26
19.532
4.3045·105
3.4465·103
TS
8
3
11.5360
5.2699·106
3.5922·104
N
x: 250.000N/m Ny: 1.000.000N/m
VNS
8
471
8.116
8.2968·106
5.4898·106
MA
8
500
8
8.2968·106
8.2815·106
AG
10
18
17.1267
7.6349·105
6.3974·104
SA
10
24
21.12
6.7663·105
2.5123·103
TS
10
3
12.4680
5.0663·105
5.1883·103
N
x: 250.000N/m Ny: 1.250.000N/m
VNS
8
118
9.528
6.6895·106
1.866·106
MA
8
500
8
6.6895·106
6.0696·106
AG
10
2
20.52
1.0673·106
1.9896·104
SA
12
3
23.672
4.5062·104
3.2251·102
TS
10
1
13.6400
5.9026·105
3.0134·103
N
x: 250.000N/m Ny: 1.500.000N/m
VNS
10
496
10.016
9.0479·105
6.8493·105
MA
10
500
10
9.0479·105
8.9861·105
AG
12
14
24.26
2.8056·105
9.3622·103
SA
10
22
26.24
5.6722·105
2.4689·103
TS
12
4
14.7000
8.7860·104
1.0494·103
N
x: 500.000N/m Ny: 500.000N/m
VNS
10
400
10.4
6.37·105
5.0705·105
MA
10
481
10.08
6.37·105
6.1494·105
AG
10
3
15.6467
5.6117·105
6.2101·104
SA
12
4
21.292
1.5327·105
1.1971·103
TS
10
3
13.9240
5.6117·105
1.0223·104
N
x: 500.000N/m Ny: 750.000N/m
VNS
12
260
12.27
8.9005·104
7.3328·104
MA
12
496
12.02
9.0984·104
8.7607·104
AG
12
4
18.9467
8.5403·104
8.4764·103
SA
16
3
26.456
4.0543·103
30.7819
TS
12
1
15.9200
1.7488·104
6.1199·102
N
x: 500.000N/m Ny: 1.000.000N/m
VNS
12
15
13.96
3.6897·104
1.3575·104
MA
12
98
13.61
3.6897·104
1.8147·104
AG
14
1
23.1267
1.1361·104
1.146·103
SA
16
3
31.672
5.6763·103
22.3771
TS
14
3
17.5640
5.2902·103
81.8767
N
x: 500.000N/m Ny: 1.250.000N/m
VNS
12
72
14.28
2.9117·104
8.6122·103
MA
12
339
12.81
2.9117·104
2.173·104
AG
16
2
27.4133
2.6399·103
2.2388·102
SA
18
2
35.388
1.23·103
4.4746
TS
16
1
19.2560
9.2541·102
23.4648
N
x: 500.000N/m Ny: 1.500.000N/m
VNS
14
363
14.636
1.1324·104
5.3271·103
MA
14
494
14.02
1.1324·104
1.0644·104
AG
18
10
30.1067
8.7844·102
85.4087
SA
16
2
38.476
9.3088·102
4.1987
TS
16
1
20.3560
8.7751·102
18.3449
N
x: 750.000N/m Ny: 750.000N/m
VNS
16
500
16
3.522·103
3.3743·103
MA
16
500
16
3.522·103
3.522·103
AG
16
15
22.7667
3.522·103
7.5575·102
SA
16
2
29.65
3.1351·103
18.3958
TS
16
16
18.9880
3.5220·103
1.4030·102
N
x: 750.000N/m Ny: 1.000.000N/m
VNS
18
500
18
1.0035·103
9.2486·102
MA
18
500
18
1.0073·103
9.8167·102
AG
18
16
27.4133
9.7849·102
1.8292·102
SA
20
4
35.736
3.3870·102
2.003
TS
18
9
21.4320
8.8434·102
24.9275
N
x: 750.000N/m Ny: 1.250.000N/m
VNS
18
38
19.848
5.2267·102
3.1132·102
MA
16
7
18.49
5.2345·102
4.3804·102
AG
20
10
30.52
3.2510·102
54.4953
SA
20
2
40.38
1.8789·102
1.3868
TS
20
2
23.7400
1.9273·102
2.8100
N
x: 750.000N/m Ny: 1.500.000N/m
VNS
18
11
20.708
3.3956·102
1.4176·102
MA
18
63
19.776
3.3956·102
1.9003·102
AG
22
16
34.3
1.0962·102
16.2374
SA
26
1
45.928
5.8144
0.032
TS
22
3
25.8480
43.3493
0.7962
N
x: 1.000.000N/m Ny: 1.000.000N/m
VNS
20
500
20
3.2888·102
3.1014·102
MA
20
500
20
3.2888·102
3.2887·102
AG
20
5
30.44
2.9623·102
44.4987
SA
20
2
38.96
2.2611·102
1.8568
TS
20
1
24.3960
1.9686·102
5.8882
N
x: 1.000.000N/m Ny: 1.250.000N/m
VNS
22
500
22
1.1748·102
1.0838·102
MA
22
500
22
1.1785·102
1.1536·102
AG
22
9
34.2467
1.1539·102
12.4818
SA
26
3
44.54
22.4217
0.1669
TS
22
3
26.7880
1.0333·102
1.7983
N
x: 1.000.000N/m Ny: 1.500.000N/m
VNS
24
500
24
46.2995
42.4606
MA
24
500
24
46.4240
45.1922
AG
24
8
37.9467
44.3415
4.5575
SA
28
7
49.4
13.0737
0.165
TS
24
1
29.6360
19.4564
0.2810
N
x: 1.250.000N/m Ny: 1.250.000N/m
VNS
24
469
24.124
46.643
41.4749
MA
24
500
24
46.643
46.6245
AG
24
3
37.0533
42.1974
4.8935
SA
26
1
46.908
13.6039
0.0983
TS
26
10
29.3760
15.1788
0.7773
N
x: 1.250.000N/m Ny: 1.500.000N/m
VNS
26
130
27.496
9.8711
7.8726
MA
26
231
27.144
10.1457
8.3843
AG
28
16
40.2333
7.8922
1.6101
SA
30
2
52.46
4.6966
0.0335
TS
28
7
32.5520
7.0430
0.2389
N
x: 1.500.000N/m Ny: 1.500.000N/m
VNS
30
500
30
3.8763
3.5545
MA
30
500
30
3.8763
3.8245
AG
30
10
42.86
3.7687
0.6741
SA
30
1
54.748
3.1197
0.015
TS
38
8
35.3000
2.3280
0.1268


### 4.1. Statistical Significance Test

In order to explore the result of simulations a statistical significance test is applied. The choice of test depends upon what it is intended to study, as Demšar [70] showed that nonparametric tests are safer and more appropriate than parametric tests for comparisons between two or more algorithms on multiple data sets.A null or no-effect hypothesis is to be formulated prior to the application of the test. It often supports the equality or absence of differences among the results of the algorithms and enables alternative hypotheses to be raised that support the opposite [71]. The null hypothesis can be represented by H0 and the alternative hypothesis by Ha. The application of the tests leads to the computation of a statistic, which can be used to reject the null hypothesis at a given level of significance α. It is also possible to compute the smallest level of significance that results in the rejection of the null hypothesis. This level is the p value, which is the probability of obtaining a result at least as extreme as the one that was actually observed, assuming that the null hypothesis is true. The use of p values is often preferred over using only fixed α levels since they provide cleaner measures of how significant the result is (the smaller the p-value, the stronger the evidence against the null hypothesis is) [72].The Friedman test [73, 74] is a nonparametric test with multiple comparisons that aims to detect significant differences between the behaviour of two or more algorithms. The null hypothesis for Friedman’s test states equality of medians between the populations. The alternative hypothesis is defined as the negation of the null hypothesis, so it is nondirectional. The first step in calculating the test statistic is to convert the original results to ranks. They are computed using the following procedure (n is the number of problems considered where i is its associated index and k is the number of algorithms included in the comparison where j is its associated index):(1)
Gather observed results for each algorithm/problem pair.(2)
For each problemi, rank values from 1 (best result) to k (worst result). Denote these ranks as rij (1≤j≤k).(3)
For each algorithmj, average the ranks obtained in all problems to obtain the final rank Rj=1/n∑irij.Thus, it ranks the algorithms for each problem separately; the best performing algorithm should have the rank of 1, the second best rank 2, and so forth. Under the null hypothesis, which states that all the algorithms behave similarly (therefore their ranksRj should be equal), the Friedman statistic Ff can be computed as(4)Ff=12nkk+1∑jRj2-kk+124which is distributed according to a chi-square distribution with (k-1) degrees of freedom.Iman and Davenport [75] derived a less conservative alternative with statistic distributed according to following F-distribution with (k-1) and (k-1)(n-1) degrees of freedom: (5)FID=n-1χF2nk-1-χF2.Table8 depicts the rank computed through the Friedman test for all heuristics considered in the exhaustive analysis (in our case of study n=21 and k=5) according to the average fitness, F¯, of each distributed Nx and Ny loading. As can be deduced from Table 8, MA with a rank of 1 is the best performing algorithm, whereas SA with a rank of 4.948 is the worst.Table 8
Average rankings of algorithms achieved by Friedman test.
Algorithm
Average rankings
MA
1
VNS
2
GA
3
TS
4.095
SA
4.948The Friedman statistic computed by (4) (distributed according to chi-square with 4 degrees of freedom) is 82.552381 and its p value is 8.256673e-11 and the Iman-Davenport extension computed by (5) (distributed according to F-distribution with 4 and 80 degrees of freedom) is 1140.526316 and its p value is 1.149264e-69. In both cases, the p value is less than the significance level α=0.05; and H0 for which there is no difference in rankings for these 5 heuristics is rejected.Rejection ofH0 must be followed by a post hoc procedure to characterize the differences between algorithms. The aim of the application of post hoc procedures is to perform a comparison considering a control method and a set of algorithms. A family of hypotheses can be defined, all related to the control method. Then, the application of a post hoc test can lead to obtaining a p value which determines the degree of rejection of each hypothesis [76]. A family of hypotheses is a set of logically interrelated hypotheses of comparisons which, in (1×n) comparisons, compares the (k-1) algorithms of the study (excluding the control) with the control method, whereas, in (n×n) comparisons, it considers the (k(k-1)/2) possible comparisons among algorithms. Therefore, the family will be composed of (k-1) or (k(k-1)/2) hypotheses, respectively, which can be ordered by its p value, from lowest to highest.In our comparison, the appropriated post hoc procedure is the Nemenyi test [77] because multiple algorithms on multiple data sets are compared. In our case of study (n=21 and k=5) there are 10 hypotheses for (5×4)/2 possible comparisons among algorithms. The p values for the Nemenyi test are shown in Table 9 ordered by p values, where in each row the first heuristic has better (lower) average ranking than the second. If the corresponding p value is smaller than α=0.05, then the first heuristic is significantly better than the second (H0 hypothesis rejected).Table 9
p values for Nemenyi test.
Hypothesis
p
N
e
m
e
n
y
i
MA vs. SA
0
MA vs. TS
0
VNS vs. SA
0
VNS vs. TS
0.000018
MA vs. GA
0.000042
GA vs. SA
0.000095
GA vs. TS
0.024796
MA vs. VNS
0.040424
VNS vs. GA
0.040424
SA vs. TS
0.097110As can be deduced from Table9, the MA proposed model significantly performs better than the others because the H0 hypothesis is rejected in all comparative pairs.


### 4.2. Analysis of Results

As indicated in Section4.1, the MA proposed model is significantly better than the others heuristics respect to the fitness function. Moreover, in many cases, the proposed model obtains a material with the lowest number of layers and such occurrences are greater than the solutions obtained by the others heuristics (see Table 7). Figure 5 shows the average number of laminas corresponding to the laminates obtained with all heuristics.Figure 5
Average number of layers for the different algorithms.As an example of results, in Table10 are shown the specifications of the best laminate obtained with the 5 heuristics in the loading case Nx=1.000.000N/m, Ny=1.500.000N/m.Table 10
Characteristics of the best laminate obtained with each algorithm.
Algorithm
Fiber
Matrix
V
f
Fitness
Stacking sequence
Thickness (mm)
MA
P-100
Peek
0.40
46.424
0
°
70
°
-
40
°
40
°
0
°
80
°
-
80
°
-
70
°
-
80
°
40
°
-
40
°
80
°
S
4.32
VNS
P-100
Peek
0.40
46.229
-
40
°
-
60
°
70
°
0
°
80
°
0
°
60
°
-
80
°
-
80
°
-
40
°
50
°
40
°
S
4.32
GA
P-100
Peek
0.40
44.341
-
50
°
-
60
°
70
°
-
10
°
80
°
0
°
60
°
-
80
°
-
80
°
-
40
°
40
°
50
°
S
4.32
SA
P-100
Peek
0.30
13.073
-
50
°
40
°
40
°
60
°
-
40
°
50
°
-
40
°
-
40
°
-
40
°
60
°
50
°
80
°
80
°
-
70
°
S
5.04
TA
P-100
Peek
0.40
19.456
50
°
-
40
°
-
50
°
-
40
°
50
°
50
°
60
°
-
60
°
-
60
°
50
°
50
°
-
50
°
S
4.32Finally, Figure6 shows graphically the different coefficients (Pki or Tsai-Wu, P1, P2, and P12) corresponding to the laminas of the best laminate obtained with the MA-based model in the loading case shown in Table 10. In this case, a large degree of uniformity of the coefficients Pki around 1 can be observed. It indicates an excellent exploitation of all the laminas.Figure 6
CoefficientsPki (Tsai-Wu), P1, P2, and P12 of each laminate in the best laminate obtained by the MA-based model.


## 4.1. Statistical Significance Test

In order to explore the result of simulations a statistical significance test is applied. The choice of test depends upon what it is intended to study, as Demšar [70] showed that nonparametric tests are safer and more appropriate than parametric tests for comparisons between two or more algorithms on multiple data sets.A null or no-effect hypothesis is to be formulated prior to the application of the test. It often supports the equality or absence of differences among the results of the algorithms and enables alternative hypotheses to be raised that support the opposite [71]. The null hypothesis can be represented by H0 and the alternative hypothesis by Ha. The application of the tests leads to the computation of a statistic, which can be used to reject the null hypothesis at a given level of significance α. It is also possible to compute the smallest level of significance that results in the rejection of the null hypothesis. This level is the p value, which is the probability of obtaining a result at least as extreme as the one that was actually observed, assuming that the null hypothesis is true. The use of p values is often preferred over using only fixed α levels since they provide cleaner measures of how significant the result is (the smaller the p-value, the stronger the evidence against the null hypothesis is) [72].The Friedman test [73, 74] is a nonparametric test with multiple comparisons that aims to detect significant differences between the behaviour of two or more algorithms. The null hypothesis for Friedman’s test states equality of medians between the populations. The alternative hypothesis is defined as the negation of the null hypothesis, so it is nondirectional. The first step in calculating the test statistic is to convert the original results to ranks. They are computed using the following procedure (n is the number of problems considered where i is its associated index and k is the number of algorithms included in the comparison where j is its associated index):(1)
Gather observed results for each algorithm/problem pair.(2)
For each problemi, rank values from 1 (best result) to k (worst result). Denote these ranks as rij (1≤j≤k).(3)
For each algorithmj, average the ranks obtained in all problems to obtain the final rank Rj=1/n∑irij.Thus, it ranks the algorithms for each problem separately; the best performing algorithm should have the rank of 1, the second best rank 2, and so forth. Under the null hypothesis, which states that all the algorithms behave similarly (therefore their ranksRj should be equal), the Friedman statistic Ff can be computed as(4)Ff=12nkk+1∑jRj2-kk+124which is distributed according to a chi-square distribution with (k-1) degrees of freedom.Iman and Davenport [75] derived a less conservative alternative with statistic distributed according to following F-distribution with (k-1) and (k-1)(n-1) degrees of freedom: (5)FID=n-1χF2nk-1-χF2.Table8 depicts the rank computed through the Friedman test for all heuristics considered in the exhaustive analysis (in our case of study n=21 and k=5) according to the average fitness, F¯, of each distributed Nx and Ny loading. As can be deduced from Table 8, MA with a rank of 1 is the best performing algorithm, whereas SA with a rank of 4.948 is the worst.Table 8
Average rankings of algorithms achieved by Friedman test.
Algorithm
Average rankings
MA
1
VNS
2
GA
3
TS
4.095
SA
4.948The Friedman statistic computed by (4) (distributed according to chi-square with 4 degrees of freedom) is 82.552381 and its p value is 8.256673e-11 and the Iman-Davenport extension computed by (5) (distributed according to F-distribution with 4 and 80 degrees of freedom) is 1140.526316 and its p value is 1.149264e-69. In both cases, the p value is less than the significance level α=0.05; and H0 for which there is no difference in rankings for these 5 heuristics is rejected.Rejection ofH0 must be followed by a post hoc procedure to characterize the differences between algorithms. The aim of the application of post hoc procedures is to perform a comparison considering a control method and a set of algorithms. A family of hypotheses can be defined, all related to the control method. Then, the application of a post hoc test can lead to obtaining a p value which determines the degree of rejection of each hypothesis [76]. A family of hypotheses is a set of logically interrelated hypotheses of comparisons which, in (1×n) comparisons, compares the (k-1) algorithms of the study (excluding the control) with the control method, whereas, in (n×n) comparisons, it considers the (k(k-1)/2) possible comparisons among algorithms. Therefore, the family will be composed of (k-1) or (k(k-1)/2) hypotheses, respectively, which can be ordered by its p value, from lowest to highest.In our comparison, the appropriated post hoc procedure is the Nemenyi test [77] because multiple algorithms on multiple data sets are compared. In our case of study (n=21 and k=5) there are 10 hypotheses for (5×4)/2 possible comparisons among algorithms. The p values for the Nemenyi test are shown in Table 9 ordered by p values, where in each row the first heuristic has better (lower) average ranking than the second. If the corresponding p value is smaller than α=0.05, then the first heuristic is significantly better than the second (H0 hypothesis rejected).Table 9
p values for Nemenyi test.
Hypothesis
p
N
e
m
e
n
y
i
MA vs. SA
0
MA vs. TS
0
VNS vs. SA
0
VNS vs. TS
0.000018
MA vs. GA
0.000042
GA vs. SA
0.000095
GA vs. TS
0.024796
MA vs. VNS
0.040424
VNS vs. GA
0.040424
SA vs. TS
0.097110As can be deduced from Table9, the MA proposed model significantly performs better than the others because the H0 hypothesis is rejected in all comparative pairs.


## 4.2. Analysis of Results

As indicated in Section4.1, the MA proposed model is significantly better than the others heuristics respect to the fitness function. Moreover, in many cases, the proposed model obtains a material with the lowest number of layers and such occurrences are greater than the solutions obtained by the others heuristics (see Table 7). Figure 5 shows the average number of laminas corresponding to the laminates obtained with all heuristics.Figure 5
Average number of layers for the different algorithms.As an example of results, in Table10 are shown the specifications of the best laminate obtained with the 5 heuristics in the loading case Nx=1.000.000N/m, Ny=1.500.000N/m.Table 10
Characteristics of the best laminate obtained with each algorithm.
Algorithm
Fiber
Matrix
V
f
Fitness
Stacking sequence
Thickness (mm)
MA
P-100
Peek
0.40
46.424
0
°
70
°
-
40
°
40
°
0
°
80
°
-
80
°
-
70
°
-
80
°
40
°
-
40
°
80
°
S
4.32
VNS
P-100
Peek
0.40
46.229
-
40
°
-
60
°
70
°
0
°
80
°
0
°
60
°
-
80
°
-
80
°
-
40
°
50
°
40
°
S
4.32
GA
P-100
Peek
0.40
44.341
-
50
°
-
60
°
70
°
-
10
°
80
°
0
°
60
°
-
80
°
-
80
°
-
40
°
40
°
50
°
S
4.32
SA
P-100
Peek
0.30
13.073
-
50
°
40
°
40
°
60
°
-
40
°
50
°
-
40
°
-
40
°
-
40
°
60
°
50
°
80
°
80
°
-
70
°
S
5.04
TA
P-100
Peek
0.40
19.456
50
°
-
40
°
-
50
°
-
40
°
50
°
50
°
60
°
-
60
°
-
60
°
50
°
50
°
-
50
°
S
4.32Finally, Figure6 shows graphically the different coefficients (Pki or Tsai-Wu, P1, P2, and P12) corresponding to the laminas of the best laminate obtained with the MA-based model in the loading case shown in Table 10. In this case, a large degree of uniformity of the coefficients Pki around 1 can be observed. It indicates an excellent exploitation of all the laminas.Figure 6
CoefficientsPki (Tsai-Wu), P1, P2, and P12 of each laminate in the best laminate obtained by the MA-based model.


## 5. Application to Specific Cases

Previously, the proposed model has been applied to general cases. In this section, the proposed model is applied to two significant specific cases comparing the results with another author’s solution.


### 5.1. Minimum Thickness Design

Le Riche and Haftka [26] developed a GA for minimum thickness composite laminate design. They consider a graphite-epoxy plate with the load case with Nx=12.500lb/in=2.189.956N/m and Ny=3.125lb/in=547.890N/m. The characteristics of the best laminated obtained by authors are shown in Table 11.Table 11
Characteristics of the best laminate obtained by Le Riche and Haftka’s GA.
Characteristic
Value
Stacking sequence
[±452/902/±453/02/±45/04/±45/02]s
N° of Laminas
48
Laminate thickness
0.240 in = 6.096 mm
F
F
m
a
x
38.002The laminate provided by the MA proposed model, under the same load case, is shown in Table12.Table 12
Characteristics of the best laminate obtained by the MA proposed model.
Characteristic
Value
Stacking sequence
[
30
4
/
-
20
2
/
-
30
4
/
20
2
/
30
2
/
-
30
2
]
s
Number of laminas
32
Laminate thickness
0.234 in = 5.94 mm
F
F
m
a
x
42.001


### 5.2. Stacking Sequence Design

Liu et al. [54] developed a permutation GA for stacking sequence design of composite laminates. They consider a graphite-epoxy laminate with the load case with Nx=-10.000lb/in=-1.751.965N/m and Ny=-2.000lb/in=-350.393N/m; Nxy=1.000lb/in=175.197N/m. The characteristics of the best laminated obtained by authors are shown in Table 13.Table 13
Characteristics of the best laminate obtained by Liu et al. permutation GA.
Characteristic
Value
Stacking Sequence
[
0
7
/
45
15
/
90
7
]
s
N° of Laminas
58
Laminate Thickness
0.290 in = 7.366 mm
F
F
m
a
x
29.002The laminate provided by the MA proposed model, under the same load case, is shown in Table14.Table 14
Characteristics of the best laminate obtained with the MA proposed model.
Characteristic
Value
Stacking Sequence
[
0
4
/
-
40
3
/
10
3
/
0
4
/
-
40
3
/
0
]
s
N° of Laminas
42
Laminate Thickness
0.275 in = 6.985 mm
F
F
m
a
x
40.010The laminate represented in Tables12 and 14 has been validated using a finite elements simulation software. To do so we employed the ANSYS 17.0 software package; additionally, the SHELL281 was used as the element of validation code, maintaining the same conditions and loadings when possible, and its characteristics can be consulted in [78].


## 5.1. Minimum Thickness Design

Le Riche and Haftka [26] developed a GA for minimum thickness composite laminate design. They consider a graphite-epoxy plate with the load case with Nx=12.500lb/in=2.189.956N/m and Ny=3.125lb/in=547.890N/m. The characteristics of the best laminated obtained by authors are shown in Table 11.Table 11
Characteristics of the best laminate obtained by Le Riche and Haftka’s GA.
Characteristic
Value
Stacking sequence
[±452/902/±453/02/±45/04/±45/02]s
N° of Laminas
48
Laminate thickness
0.240 in = 6.096 mm
F
F
m
a
x
38.002The laminate provided by the MA proposed model, under the same load case, is shown in Table12.Table 12
Characteristics of the best laminate obtained by the MA proposed model.
Characteristic
Value
Stacking sequence
[
30
4
/
-
20
2
/
-
30
4
/
20
2
/
30
2
/
-
30
2
]
s
Number of laminas
32
Laminate thickness
0.234 in = 5.94 mm
F
F
m
a
x
42.001


## 5.2. Stacking Sequence Design

Liu et al. [54] developed a permutation GA for stacking sequence design of composite laminates. They consider a graphite-epoxy laminate with the load case with Nx=-10.000lb/in=-1.751.965N/m and Ny=-2.000lb/in=-350.393N/m; Nxy=1.000lb/in=175.197N/m. The characteristics of the best laminated obtained by authors are shown in Table 13.Table 13
Characteristics of the best laminate obtained by Liu et al. permutation GA.
Characteristic
Value
Stacking Sequence
[
0
7
/
45
15
/
90
7
]
s
N° of Laminas
58
Laminate Thickness
0.290 in = 7.366 mm
F
F
m
a
x
29.002The laminate provided by the MA proposed model, under the same load case, is shown in Table14.Table 14
Characteristics of the best laminate obtained with the MA proposed model.
Characteristic
Value
Stacking Sequence
[
0
4
/
-
40
3
/
10
3
/
0
4
/
-
40
3
/
0
]
s
N° of Laminas
42
Laminate Thickness
0.275 in = 6.985 mm
F
F
m
a
x
40.010The laminate represented in Tables12 and 14 has been validated using a finite elements simulation software. To do so we employed the ANSYS 17.0 software package; additionally, the SHELL281 was used as the element of validation code, maintaining the same conditions and loadings when possible, and its characteristics can be consulted in [78].


## 6. Conclusions

In this paper, we have shown that the use of techniques based on systematic exploitation of knowledge of the problem, together with the combination of metaheuristics based on population as well as in local search, offers more effective and dynamic methods for the solution of laminated polymers design. To this end, a study of the main techniques used in the design of laminated polymers was firstly performed, in order to carry out a comparative analysis.A Memetic Computing-based model has been presented for the design of Symmetric Laminated Composites and Structures. This model implements a general encoding for the design of composites and a fitness function that has taken into account economic and safety criteria in design. Also, a dynamic reproductive operator is presented, in which the classic crossover operator is modified in order to improve the dynamics in the change of the number of layers.Finally, a set of local search operators are implemented. These allow a lamina to be added or removed, to change its orientation, type of matrix, type of fiber, or volume fraction. Two new operators have been defined and added to the set to improve the exploitation of the solutions in the region of the search space.The proposed Memetic Computing model has been subjected to a broad analysis and applied in two specific cases. Firstly, the model has been applied to the design of a plate under distributedNx and Ny loading, and the results have been compared with those obtained by other well-known design methods. In most cases, the minimum number and average number of laminas are lesser. Also, the proposed model has been compared with two significant models found in the literature, obtaining comparable results. These results show that the proposed model is a general and flexible design method for Symmetric Laminated Composites and Structures. We consider that these results can be applied to real-life scenarios with a high reliability. Finally, all solutions obtained by the MA proposed model in Sections 5.1 and 5.2 have been validated by the ANSYS software package using the same conditions and load system.



---

*Source: 4723863-2017-09-28.xml*