Multi-Objective High-Fidelity Optimization Using NSGA-III and MO-RPSOLC

N.Ganesh,Uvaraja Ragavendran,Kanak Kalita,Paras Jain and Xiao-Zhi Gao

1Department of Computer Science and Engineering, Vel Tech Multi Tech Dr.Rangarajan Dr.Sakunthala Engineering College,Chennai,600062,India

2Department of Electronics and Telecommunication Engineering, MPSTME SVKM’S Narsee Monjee Institute of Management Studies,Shirpur,425405,India

3Department of Mechanical Engineering,Vel Tech Rangarajan Dr.Sagunthala R&D Institute of Science and Technology,Avadi,600062,India

4School of Computing Science and Engineering,VIT Bhopal University,Sehore,466114,India

5School of Computing,University of Eastern Finland,Kuopio,70211,Finland

ABSTRACT Optimizing the performance of composite structures is a real-world application with significant benefits.In this paper, a high-fidelity finite element method (FEM)is combined with the iterative improvement capability of metaheuristic optimization algorithms to obtain optimized composite plates.The FEM module comprises of ninenode isoparametric plate bending element in conjunction with the first-order shear deformation theory(FSDT).A recently proposed memetic version of particle swarm optimization called RPSOLC is modified in the current research to carry out multi-objective Pareto optimization.The performance of the MO-RPSOLC is found to be comparable with the NSGA-III.This work successfully highlights the use of FEM-MO-RPSOLC in obtaining highfidelity Pareto solutions considering simultaneous maximization of the fundamental frequency and frequency separation in laminated composites by optimizing the stacking sequence.

KEYWORDS Composites; finite element; optimization; Pareto; swarm intelligence

1 Introduction

Composite structures, owing to their superior properties like excellent strength to weight ratios and superb stiffness to weight ratios are steadily replacing steel and other metallic structures in structural engineering.As such, these composite structures often operate in various static and dynamic loading conditions [1].Often these composite structures carry vibrating machinery/equipment on them.To avoid resonance, it becomes necessary that these structures are designed in such a way that the vibrating equipment operates well outside the first few prominent natural frequencies of the structure.

Simulating the static or dynamic behaviour of laminated composites using numerical methods like finite element analysis (FEA)has achieved a lot of attention.Such high-fidelity numerical methods can be coupled with optimization methods to predict the optimal parameter combinations.For a given design or a specific application, often the geometric parameters of composite plates are unalterable.In such cases, altering the stacking sequence of laminated plates to optimize its performance is the most convenient approach.However, owing to the wide design space, i.e.,ply angle range of ±90?for each lamina, it is physically impossible to experimentally obtain the optimal combination of ply angles, especially when the number of plies is high.For example, for a symmetric laminated composite plate having 4 layers, i.e., 2 design variables, the total number of possible ply angle combinations is 1812(considering discrete ply angles with 1?increment).Thus,it is clear that the complexity of the search space and its dimension will increase substantially if the number of layers in the composite plate is increased.Adali et al.[2] used a mathematical programming approach to predict the optimal fundamental frequencies in symmetric laminates.Narita [3] used a Ritz-based layerwise optimization approach, which was able to achieve appreciable results with low computational costs.However, as pointed out by Apalak et al.[4], such mathematical methods are often sensitive to gradients and sometimes may get trapped in local optima.Diaconu et al.[5] in a study on thick laminated plates obtained the optimal lamination parameters and corresponding ply angles without encountering any multi-modality.

To address the shortcomings of gradient-based optimizers, many researchers have relied on metaheuristic optimization techniques.Le Riche and Haftka were among the pioneering researchers to use a genetic algorithm (GA)in laminate optimization.They carried out buckling load maximization [6], thickness minimization [7], stacking optimization [8], etc.Apalak et al.[4]carried out stacking sequence optimization to maximize fundamental frequency using an artificial bee colony (ABC)and GA.Hemmatian et al.[9] used an imperialist competitive algorithm (ICA)along with GA and ant colony optimization (ACO)to simultaneously optimize the weight and cost of a rectangular composite plate.They reported that ICA outperformed GA and ACO.Kalita et al.[10] used GA, PSO (particle swarm optimization)and Cuckoo search to simultaneously optimize fundamental frequency and separation between the first two frequencies.However, both Hemmatian et al.[9] and Kalita et al.[10] used a weighted-sum approach to multi-objective optimization.Correia et al.[11] considered stacking sequence as the design variables while maximizing frequency parameters and minimizing strain energy in composite plates with piezoelectric layers.Vo-Duy et al.[12] used non-dominated sorting genetic algorithm II (NSGA-II)in conjunction with a finite element analysis to minimize weight and maximize the frequency of composite plates.Ghasemi et al.[13] used a similar NSGA-II based strategy to minimize the cost of composite shells while increasing their buckling strength.An et al.[14] studied multi-objective optimal designs of hybrid laminates using genetic algorithms.Serhat et al.[15] carried out fundamental frequency, buckling load and effective stiffness maximization using lamination parameters.Serhat et al.[16] also carried out multi-objective optimization of stiffened composite fuselages.

The literature search reveals that though a lot of work on metaheuristic and FEA-based high-fidelity single-objective optimization has been carried out, similar high-fidelity optimization in a multi-objective perspective is rare.A novel particle swarm optimization called RPSOLC(repulsive particle swarm optimization with local search and chaotic perturbation), demonstrated by the author(s)as an effective single optimization tool in their previous works [10,17,18] has been further enhanced to carry out multi-objective Pareto optimization.By comparing with a powerful and recent multi-objective algorithm called the NSGA-III (non-dominated sorting genetic algorithm-III), the efficacy and utility of MO-RPSOSLC in designing laminated composites are demonstrated.

2 Methodology

The methodology used in the current research work is shown in Fig.1.Laminated composite plate design considering multiple natural frequency criteria is considered as the design problem.The design problem is introduced in detail in Section 2.1.The ply orientations of the individual lamina are selected as the design variables.A finite element module is developed in Fortran language that can accurately determine the natural frequencies of the laminated plate.The finite element formulation is discussed in detail in Section 2.2.Next, two independent multi-objective optimization modules (NSGA-III and MO-RPSOLC)for carrying out Pareto optimization are developed in Fortran.These are discussed in detail in Sections 2.3 and 2.4, respectively.The FEM module is incorporated as the objective function in the NSGA-III and MO-RPPSOLC modules.Over several generations, the Pareto solutions are iteratively improved by the optimization techniques.The Pareto fronts by the NSGA-III and MO-RPSOLC are compared and the superior one is selected.Further, suitable compromise solutions relating to various test cases are selected by using a multi-criteria decision-making method called TOPSIS.The details of the TOPSIS are discussed in Section 2.5.

Figure 1: The methodology used in the current research work

2.1 Problem Statement

In this paper, laminated composite plates of the configuration as shown in Fig.2 is used.An eight-ply squaresymmetric Graphite-epoxy [19] composite laminate having a thickness(h=0.01)and material properties as follows is considered:

E1=138GPa,E2=8.96GPa,G12=7.1GPa,υ12=0.3.

The boundary condition at the edges of the plates may be clamped (C), simply supported(S)or free (F).Based on these conditions, 18 different boundary condition combinations are considered for the laminated plates.

Figure 2: Typical laminated plate considered in the study

The fundamental frequency(λ1)and the difference between the first two natural modes(λ21)are to be maximized simultaneously.The stacking sequence (i.e., the ply angles denoted byθ)of the laminated composites are considered as the design variables to achieve this multi-objective optimization.Optimizing the stacking sequence of laminated composite plates to optimize its performance is a combinatorial optimization problem, which may be stated as,

As such two powerful metaheuristic optimization algorithms—NSGA-III and MO-RPSOLC are used.Both these algorithms independently generate a set of non-dominated solution known as the Pareto front.Each solution within the Pareto front is a compromise multi-objective solution.Since the frequency parameters are quite sensitive to the stacking sequence, a high-fidelity finite element analysis is carried out which runs in conjunction with the optimization algorithms.

2.2 Finite Element Method Formulation

In this work, high-fidelity structural analysis of the laminated composite plates is carried out using a highly accurate finite element analysis (FEA)module.In past, this FEA module has been used by Kalita et al.[1] to accurately compute the Eigen frequencies in taper plates, laminated plates [20] and isotropic plate with holes [21].

The FEA module is based on a nine-node isoparametric plate bending element.The element contains four corner nodes, four mid-edge nodes and one centre node.The effect of shear deformation is an important consideration for laminated composites since they are weak in shear.This is included in the current FEA module by assuming the considerations of FSDT (First-order Shear Deformation Theory)which states that the normal to the midplane would remain straight after deformation but may or may not remain normal.

whereφxandφyrepresent the average shear rotation over the entire plate thickness andθx, whileθyare denote total rotations in bending.u,vare the in-plane displacements whereaswis the transverse displacement.

The nine-node isoparametric plate bending element maps the laminated plate by using the following equations:

where(x,y)is the location of any given point within the element, therth nodal point’s location within the element is given by(xr,yr)and the assumed Lagrangian interpolation function is denoted asNr.

For any composite laminate, the stress-strain relation can be expressed as,

{σ} is the generalized stress vector.

where the in-plane force resultants are represented asNx,Ny,Nxy; bending moments in thexandydirections are indicated asMx,My; the twisting moment resultant is indicated asMxy.Qx,Qyare the transverse shear force resultants.

In terms of displacement fields, the strain vector is,

Aij,BijandDijrepresent the extensional stiffness, extensional-bending stiffness and bending stiffness coefficients respectively.The number of individual laminas forming the composite plate is indicated asn.The material coefficients of the kth layer are denoted as

whereE1andE2are the longitudinal and transverse Young’s modulus, respectively.μ12andG12,G13,G23are the in-plane Poisson’s ratio and shear moduli, respectively.

The strain vector can be summarized as,

where [B] is the strain matrix, {δ}is the nodal displacement vector and [K] is the stiffness matrix.|J|is the Jacobian matrix.

The lumped mass matrix is computed as,

where,

[Nu]=[[Nr][N0][N0][N0][N0]]

[Nv]=[[N0][Nr][N0][N0][N0]]

[Nw]=[[N0][N0][Nr][N0][N0]]

where [N0] is a null matrix of the order 1×9.

[Nu]T[Nu] and [Nv]T[Nv] are responsible for the in-plane movements of mass.[Nw]T[Nw] is related to inertia is the transverse movement of mass.andare responsible for rotary inertia.

[K0] and [M0] are the overall stiffness matrix and mass matrix respectively.The plates’equation of motion is,

whereωis the frequency of the plate, calculated by solving Eq.(19)using the simultaneous iterative technique.In this article, all the frequency parameters are expressed in non-dimensional form as,

Kalita et al.[1,21] in their previous works has shown that the current FEM formulation can determine the natural frequencies with less than 1% error as compared to higher-order shear deformation theory and analytical solutions for thin and moderately thick plates.A mesh size of 18 × 18 was found to be satisfactory in those studies.However, in the current work, to keep the computational time manageable, a reduced mesh of 4 × 4 is used as shown in a previous optimization study [18] by the author(s).

2.3 Optimization with NSGA-III

In this work, the non-dominated sorting genetic algorithm III (NSGA-III)[22,23] is used for carrying out Pareto optimization.The multi-objective optimization problem is stated in Eq.(1).The NSGA-III is based on Darwin’s theory of natural evolution.The algorithm starts by assuming a set (called population)of possible solutions (called individuals).Each individual has certain properties (called chromosomes)which can be combined with others, mutated and thus, improved.These improvements in the fitness of the individuals go on until a termination criterion is met.The FEA-NSGA-III used in the current research is realized by using the pseudo-code given in Algorithm 1.

Algorithm 1: FEA-NSGA-III BEGIN Objective function, λ=f (x),x=(θ1,...,θd)Initiate generation, t=0 Initiate a random population Pt of size npop by encoding Calculate fitness of Pt using FEM DO{Create a new population Ct from Pt by crossover Modify the population Ct by mutation Calculate the fitness of Ct using FEM(Continued)

Algorithm 1.(continued)Set Rt=Pt ∪Ct Apply non-dominated sorting on Rt and find F1, F2,...St={},i=1;FOR |St|≤npop DO{St=St ∪Fi i=i+1}IF |St|=npop DO Pt+1=St; BREAK ELSE{Pt+1=∪l?1 j=1Fj Normalize St using min and intercept points of each objective Associate each member of St to a reference point Choose npop ?|Pt+1|members from Fl by niche preserving operator}t=t+1}UNTIL (t=tmax)report the Pareto front END

2.4 Optimization with MO-RPSOLC

In this article, a novel memetic version of PSO (Particle Swarm Optimization)called RPSOLC(Repulsive Particle Swarm Optimization with Local Search and Chaotic Perturbation)is used for multi-objective optimization of laminated composites.In past, single-objective RPSOLC has been used by Kalita et al.[17] for optimizing the laminated composites.Kalita et al.[10] has also used RPSOLC to conduct weighted sum multi-objective optimization.However, this the first instance of its usage in the Pareto optimization scenario.The architecture used for selecting the Pareto optimal solutions in the current MO-RPSOLC (multi-objective RPSOLC)is similar to that used by Coello et al.[24].The working of the RPSOLC algorithm is explained below:

Any typical PSO algorithm [25] starts by generating a set of random solutions (called particles)in the design space.In PSO terminology, this set of solutions is collectively called the swarm.Each particle within the swarm is aware of its personal best position called thepBest.The particles in the swarm are also aware of the overall best position (or solution)found so far called thegBest.All particles try to move towards thegBestposition by using the following two rules to update their velocity and position.

wherevtandvt+1are the velocity in the current and the next generation, respectively.Similarly,xtandxt+1are the positions of the particle in the current and the next generation, respectively.r1andr2are two random number between 0 to 1.c1andc2are called the cognitive and social parameters.The effect of the velocity of the previous generation on the velocity of the current generation is controlled by the inertia weightω.

Due to the traditional PSO’s tendency to get stuck in local optima, Urfalioglu [26] proposed the RPSO in which the particles update their velocity using the relation below:

whereα,β,γare constants.The termα.r1.(pBest?xt)leads any particle towards its self-best position.The termω.β.r2.leads the particle away from a randomly chosen particle from the swarm.The termω.γ.r3.vtris responsible for the exploration of new search regions.

The traditional RPSO was further upgraded by incorporating the memetic attributes suggested by Santos et al.[27,28].Each particle is endowed with the ability to conduct a local search by visiting its surroundings.The visit is independent of the gradient in any directions and thus, free from bias.If any particle is seen to be trapped in one position even after a pre-defined number of generations, a small random disturbance in its velocity (called the chaotic perturbation)is inserted.

whererchaosis the chaotic perturbation.This helps the RPSOLC algorithm to avoid getting trapped in local optima pits.The pseudocode of the current FEA-MO-RPSOLC approach is shown in Algorithm 2.

Algorithm 2: FEA-MO-RPSOLC BEGIN Objective function, λ=f (x),x=(θ1,...,θd)Initiate generation, t=0 Initiate a random swarm SWt of npop particles Calculate the fitness of SWt using FEM Record non-dominated particles in an external repository Update pBest Update gBest DO{FOR each particle i DO{Update velocity of particles using Eq.(23)Update position of particles using Eq.(24)Update pBest}Update gBest Update external repository t=t+1}(Continued)

Algorithm 2.(continued)UNTIL (t=tmax)report the Pareto front END

2.5 Multi-Criteria Decision Making with TOPSIS

For an MCDM problem consisting ofmalternatives andncriteria, letD=xijbe a decision matrix, wherexij∈R.The Pareto fronts obtained from NSGA-III and MO-RPSOLC are used as the decision matrix in this paper.

There are different objective and subjective weight allocation methods.In this paper, the standard deviation approach described in the preceding section is used.The weight vector may be expressed as,

A normalized decision matrixnijof each criterion is created using Eq.(27),

Next, the weighted normalized matrix is created using Eq.(28),

The ideal positive (best)and ideal negative (worst)solutions are then estimated using Eqs.(29)and (30), respectively.

where B is a vector of benefit function and C is the vector of the cost function, fori∈「1,m?andj∈「1,n?.

The separation measurement and relative closeness coefficient are then determined.In TOPSIS, the difference of each response from the ideal positive (best)solution is given by Eq.(31).

fori∈[1,m] andj∈[1,n]

Figure 3: Pareto front by NSGA III and MO-RPSOLC for (a)CCCC (b)CCCF (c)CCCS (d)CCFF

Similarly, the difference between each response from the ideal negative (worst)solution is given by Eq.(32).

fori∈[1,m] andj∈[1,n]

Figure 4: Pareto front by NSGA III and MO-RPSOLC for (a)CFCF (b)CFFF (c)CSCF (d)SCCF

The corresponding closeness coefficient (CCi)of theithalternative is calculated using Eq.(33).

where 0 ≤CCi≤1,i∈[1,m]

The final step is to rank the alternatives in decreasing order of closeness coefficient value.

Figure 5: Pareto front by NSGA III and MO-RPSOLC for (a)SCFF (b)SCSC (c)SCSF (d)SFCF

3 Results and Discussion

Figs.3–7 show the high-fidelity Pareto fronts generated by the NSGA-III and MO-RPSOLC algorithms for 18 different boundary conditions.It is seen that despite the geometric configuration and the material of the laminated composites being the same, in these 18 cases, the boundary condition plays a significant role in the shape of the Pareto fronts.For instance, in cases like CCCC, CCSS, SSSS, etc.the Pareto fronts are discontinuous whereas, in the case of boundary condition combinations like CCCF, CSCF, SCSC the Pareto fronts are smooth and continuous.The discontinuity in Pareto fronts has perhaps arisen due to the algorithm encountering the concave portion of the objectives.This has led to the presence of regions of the objective functions where they are dominated.The severe discontinuity in Pareto fronts in CCCC SSCC,SSFF and SSSS boundary conditions indicates that the presence of several distinct feasible regions with gaps where the solutions are dominated.

Figure 6: Pareto front by NSGA III and MO-RPSOLC for (a)SSCC (b)SSCF (c)SSFF (d)SSSC

Figure 7: Pareto front by NSGA III and MO-RPSOLC for (a)SSSF (b)SSSS

The optimal stacking sequences corresponding to the limits of the Pareto fronts for different boundary condition cases are reported in Tab.1.It is seen that in general, the limits of the Pareto fronts for both NSGA-III and MO-RPSOLC show negligible variation.However, it is interesting to observe that the optimal ply angles reported by NSGA-III and MO-RPSOLC mostly vary by either (approx.)5?or in the orientation direction (±)of the fibre.

The performance of the NSGA-III and MO-RPSOLC in terms of the spread of the Pareto fronts for all the 18 combinations of boundary conditions is shown in Fig.8.It is seen that in general, the NSGA-III has better uniformity in the Pareto solution spread as compared to the MO-RPSOLC.In most cases, the MO-RPSOLC solutions spread is seen to be skewed.This in contrast to the NSGA-III’s uniformity in the spread, indicates that perhaps the MO-RPSOLC has not explored the design space sufficiently.However, in all the cases the extremum points of the Pareto fronts for both the NSGA-III and MO-RPSOLC seems to be similar.It is also evident from Figs.8a, 8c and 8e that as the constraints increase the average (of the Pareto front range)fundamental frequency increases.This is because the increase in constraints (i.e.,restrictions on the degree of freedoms)causes an increase in the stiffness of the structure which,in turn, is responsible for the increased fundamental frequency.As seen from Figs.8b, 8d and 8f,baring a few exceptions like CCCC, SSCC boundary condition combinations, as the constraints increase, the range of the frequency separation between the first two natural frequencies also increase.It is worth pointing out here that the average time taken for MO-RPSOLC simulations is approximately 2/3rd of the time taken by NSGA-III.On average, it took 265 min to carry out each NSGA-III simulation on a Dell Inspiron 15-3567 series windows system with Intel(R)CoreTMi7-7500U CPU @ 2.70 GHz, Clock Speed 2.9 Ghz, L2 Cache Size 512 and 8 GB ram.

Figure 8: Spread of the Pareto fronts for all 18 boundary conditions

TOPSIS is used to select the best solutions from the Pareto fronts based on three different test cases.In test case 1, only the ranked 1 solution by TOPSIS is considered; in test case 2,the top 3 ranked solutions are considered.Similarly, for test case 3, the top 10 ranked TOPSIS solutions from the combined Pareto fronts are considered.Additionally, to draw unbiased and comprehensive conclusions based on the TOPSIS analyses, the criteria weightw1is varied from 0.01 to 0.99 with 0.01 increments for all the boundary conditions.Based on the considered test cases, the optimum predictions for eachw1is recorded.

The percentage share of the NSGA-III and MO-RPSOLC solutions in the TOPSIS selected optimal solutions (99 in total, 1 solution for eachw1considering test case 1 is calculated and shown in Fig.9a.It is seen that, in general, MO-RPSOLC solutions have a higher share in the TOPSIS selected optimal solutions.However, in the instances of SSCC, SSSC, SSCF, SSSF,SCCF and SCFF boundary conditions NSGA-III solutions find better representation.Similarly,in Fig.9b, the TOPSIS selected the top 3 optimal solutions for eachw1is recorded.The share of NSGA-III solutions in the TOPSIS selected optimal solutions is seen to have significantly improved as compared to test case 1.Similar improvement for NSGA-III is seen in Fig.9c for test case 3 considering TOPSIS selected top 10 optimal solutions.

Figure 9: The percentage share of NSGA-III and MO-RPSOLC in TOPSIS selected best solutions considering (a)test case 1: ranked 1 solution (b)test case 2: ranked 1–3 solutions (c)test case 3:ranked 1–10 solutions

4 Conclusion

In this research, two popular nature-inspired multi-objective optimization algorithms are evaluated based on their performance in a high-fidelity application.A combinatorial optimization problem of maximizing the first frequency and the difference between the first two natural frequency by optimizing the stacking sequence is chosen.The optimization problem is carried out for 18 different boundary condition combinations.The performance of the algorithms is compared based on the spread of their Pareto solutions as well as by using a novel TOPSIS-based analysis.The comprehensive analysis carried out on the 18 different test cases indicates the MO-RPSOLC to be on par with the reliable and widely used NSGA-III.It is worth mentioning that this is the first instance of using a high-fidelity finite element analysis in conjunction with NSGA-III or MO-RPSOLC to design an optimized structure.Based on the encouraging results of the article it can be concluded that this methodology can be effectively used to predict highly accurate and optimized structural designs which would lead to significant savings in terms of prototype testing.

Acknowledgement:Authors acknowledge the support of their respective institutes.

Funding Statement:The authors received no specific funding for this study.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.