Metaheuristic Optimization of Laminated Composite Plates with Cut-Outs

: The stacking sequence optimization of laminated composite plates while maximizing the structural performance or minimizing the weight is a subject investigated extensively in the literature. Meanwhile, research on the optimization of laminates with cut-outs has been relatively limited. Cut-outs being an indispensable feature of structural components, this paper concentrates on the stacking sequence optimization of composite laminates in the presence of circular cut-outs. The buckling load of a laminate is used as a metric to quantify the structural performance. Here the laminates are modeled as carbon ﬁber-reinforced composites using the ﬁnite element analysis software, ABAQUS. For the optimization, the widely used harmony search algorithm is applied. In terms of design variables, ply thickness, and ﬁber orientation angles of the plies are used as continuously changing variables. In addition to the stacking sequence, another geometric variable to consider is the aspect ratio (ratio of the length of the longer sides to the length of the shorter sides of the plate) of the rectangular laminates. The optimization is carried out for three different aspect ratios. It is shown that, by using dispersed stacking sequences instead of the commonly used 0 ◦ / ± 45 ◦ / ± 90 ◦ ﬁber angle stacks, signiﬁcantly higher buckling loads can be achieved. Furthermore, changing the cut-out geometry is found to have a signiﬁcant effect on the structural performance.


Introduction
Fiber-reinforced composites are frequently used in structural components due to their high strength-to-weight ratio and stiffness properties. Laminated composite plates, in particular, have been widely applied in civil, mechanical, aerospace, and marine engineering [1][2][3][4][5]. A novel application of composite laminates has been in the field of structural retrofitting, particularly the retrofitting of reinforced concrete and masonry structures, as presented by Mohee et al. [3,4] and Noor-E-Khuda et al. [5]. Because of the advantages of using laminated composites in structural applications, their performance has been researched extensively. In these fields, various optimization studies have been carried out to maximize the performance of these structures while reducing cost [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. Vo-Duy et al. [24,25] presented the weight minimization and natural frequency maximization problems of laminated plates using an adaptive elitist differential evolution algorithm treating the ply thicknesses, fiber orientations, and the fiber volume fractions as the variables of optimization. As the constraint of optimization, the first fundamental frequency of the laminated plate was kept above a certain lower bound in every iteration. It was observed that for convergence to a minimum weight, around a thousand iterations were necessary with the elitist differential evolution algorithm. Le-Anh et al. [26] investigated the strain energy minimization and fundamental frequency maximization of folded laminated composite plates using an adjusted differential evolution algorithm. The design variables were chosen to be the fiber orientation angles which were constrained to have integer values in the range [−90 • , 90 • ]. Keshtegar et al. [27] worked on buckling load maximization using an improved partial swarm optimization algorithm. The number of layers varied from 2 to 10. In addition to uniaxial loading, also biaxial loading case was analyzed. The present study focuses on the effects of introducing cut-outs on the structural performance of laminated composite plates.

Prior Research on Buckling Analysis of Laminated Composite Plates with Cut-Outs
Cut-outs in structural components are inevitable and serve various purposes, such as reducing the weight of the component like castellated beams or providing access to other parts of the structure [28]. In the present study, the structural performance is quantified using the buckling load of a plate under in-plane compressive forces. Buckling is a phenomenon associated with the reduction in structural stiffness. This reduction in stiffness leads to rapid growth in the displacements, which should be prevented. Especially in thin-walled structures, it is important to maximize the buckling load to stay in the pre-buckling regime throughout the lifespan of a structure. Buckling load, natural frequency, and maximum displacement under bending load are commonly used measures of structural performance. Since laminated plates made of carbon fiber composites are very thin, lightweight structures, they are prone to buckling. Therefore, buckling load is used as a primary indicator of carbon fiber-reinforced polymer (CFRP) laminates' performance. Since the sequence of ply thicknesses and fiber orientation angles of the plies have the most significant effect on the buckling load, this research aims to find the optimum stacking sequences for different plate aspect ratios. The harmony search method, a well-established and very effective metaheuristic optimization algorithm, is employed to achieve this goal. CFRP is selected as the plate material due to its widespread use in industry. Before presenting the methodology and findings of the present work, some of the notable studies in laminated composite plates with cut-outs are summarized in the following section.
Anish et al. [29] investigated the uniaxial and biaxial buckling of laminated composite plates featuring cut-outs by applying the improved shear deformation theory. It was found that the stacking sequence has a significant effect on the mode of deformation. Komur et al. [30] analyzed the buckling response of woven-glass-polyester laminated composite plates with elliptical cut-outs using the finite element method (FEM). They performed an extensive parametric study to investigate the shape variation and the position of the elliptical cut-out. After the cross-ply [(0 • /90 • )2]s and angle-ply s of stacking sequences were simulated, the cross-ply stacking sequence was found to perform better than the other stacking sequences. Barut et al. [31] presented a semi-analytical method for predicting the buckling load of a laminate with a cut-out. They investigated a square laminate with a central circular cut-out subject to uniform end-shortening. The diameter of the cut-out was chosen to be 60% of the plate side length. The thickness-to-width ratio of the laminate was found to be less than 0.05, such that elastic buckling would be expected to occur before the first-ply failure. Rajanna et al. [32] studied the effect of circular cut-outs and non-uniform in-plane edge loads on the buckling response of laminated panels. Based on finite element analyses, it was found that plates with larger cut-outs perform better than plates with smaller cut-outs when reinforced with ring stiffeners. Furthermore, they observed that panels with (±45 • )s layup performed better than other stacking sequences. Bhardwaj et al. [33] investigated the vibrational response of laminates with triangular cut-outs, focusing on the effects of geometric parameters such as aspect ratio, number of layers, fiber orientation angle, and distance between cut-outs. They observed that the natural frequency of the laminate was proportional to the number of layers and inversely proportional to the cut-out size. The fully clamped boundary condition resulted in the highest frequencies. Lahouel and Guen-foud [34] performed a vibrational and buckling analysis of laminated composite plates with and without rectangular cut-outs under uniaxial compression. It was concluded that a fiber orientation angle of 45 • is most favorable in terms of buckling load and corresponds to the highest natural frequencies. Similarly, Brethee [35] carried out vibration analysis for laminated composite plates with rectangular and circular cut-outs at the plate center and concluded that natural frequencies are not only dependent on the size and shape of the cut-out but also tended to increase with an increase in cut-out size. In another study, Mahdavi et al. [36] examined the thermal buckling load of a square hybrid composite plate with a quasi-square cut-out. A genetic algorithm was employed to maximize the critical buckling temperature by varying certain design parameters, such as the bluntness of the cut-out corners and the cut-out size-to-plate size ratio. The study concluded that plates with quasi-square cut-outs had higher thermal buckling temperatures than plates with circular cut-outs. Furthermore, it was found that the thermal buckling temperature increased with the cut-out size to plate length ratio. Plates with quasi-square cut-outs were further studied by Moussavian and Jafari [37]. Sahu and Datta [38] investigated the dynamic stability of laminated composite curved panels featuring cut-outs, finding that the effect of curvature in laminated composite curved panels is inversely proportional to the cut-out size. Narayana et al. [39] analyzed the thermal buckling of a laminated graphite/epoxy composite plate featuring rectangular and elliptical cut-outs, investigating the effect of the cut-outs on the thermal buckling behavior of 16-ply laminate for seven different stacking sequences. It was found that the magnitude of the thermal buckling temperature increased with the cut-out size-to-plate size ratio, and the orientation angle of the cut-out had a significant effect on the thermal buckling temperature. Another study involving the thermal buckling of rectangular laminated composite plates with circular cut-outs was carried out by Shaterzadeh et al. [40]. In addition to examining plates with central cut-outs, the effect of locating the cut-out at different parts of the plate (i.e., at locations other than the center of the plate) was studied. Their notable conclusions were that a central cut-out increases the thermal buckling temperature and that moving the cut-out from the center of the plate towards the clamped edges leads to additional increases in the buckling temperature. Topal and Uzman [41] studied symmetrically laminated composite plates with central circular cut-outs using the modified feasible direction optimization method. The main outcome was that an increase in the number of plies has a favorable effect on the buckling load up to a certain point, beyond which increasing the ply count further has only a minimal effect on the buckling load.

Scope and Objectives
Although the stacking sequence optimization of laminated composite plates is a subject that has received considerable attention within the literature, relatively limited research has been conducted on the optimization of laminates with cut-outs. The present study focuses on obtaining stacking sequences that maximize the buckling load of CFRP laminates with circular cut-outs using the harmony search metaheuristic algorithm. Plates of three different aspect ratios are analyzed with the finite element analysis software, ABAQUS. The design variables in the optimization process are ply thickness and fiber angle orientation. Primarily, the laminates consist of nine layers, as per Führer [42]. Concerning the choice of optimization technique, Cakiroglu et al. [23] have shown that the harmony search (HS) optimization technique was capable of obtaining dispersed stacking sequences that perform better than conventional sequences, such as the [45, −45, 0, 90, 0, 90, 0, −45, 45] fiber angle sequence used by Führer [42]. The present study adopted the optimum stacking sequence from the analyses done by Cakiroglu et al. [23] as a starting configuration and further improved this stacking sequence in the presence of cut-outs. In addition to circular cut-outs, the use of rectangular and triangular cut-outs is investigated, including an analysis of the effects of each on the structural performance. The novelty of this current study is that unlike many other studies in the literature the fiber orientation angles are allowed to take any value between 0 • and 90 • . Considering the prevalence of cut-outs in composite laminates, choosing the most efficient stacking sequence bears great importance.

Optimization Procedure
In this section, the optimization method and subsequent finite element modeling procedure will be discussed. The stacking sequences are optimized using the harmony search technique to obtain the best possible performance from plates with cut-outs. As the performance indicator of a plate, its buckling load is used. This buckling load is obtained from an eigenvalue buckling analysis procedure performed using ABAQUS. The optimization procedure is carried out for three different aspect ratios (the ratio of the longer side length to the shorter side length) of a rectangular CFRP laminate with circular cut-outs. Plates with square and triangular cut-outs are also simulated to show the effect of cut-out shapes on the structural behavior. In these simulations, the cut-outs are introduced in such a way that all of them have equal area and position of the centroid. Each one of the simulated laminates contains a certain number of cut-outs according to the length of the longer side of the plate such that the plate with the aspect ratio a/b = 3 contains 3 cut-outs whereas the plate with a/b = 1 contains one cut-out. Each one of these cutouts has a diameter equal to one-third of the shorter side length of the plate. All plate configurations have the same length for their shorter sides, which is 150 mm long, resulting in a cut-out area of 1964 mm 2 . As the starting configuration for the stacking sequence, the best-performing layup is taken from the previous study by the authors for all three aspect ratios [23]. These layup configurations are further optimized, and the results are compared to the performance of plates with cut-outs having the regular stacking sequence from Führer [42].

Harmony Search Optimization Process
In this Meta-heuristic optimization, advanced algorithms have been increasingly applied to science and engineering in recent years. The harmony search method is one of the most successful algorithms developed in the field of metaheuristic optimization. This method was invented by Geem et al. [43] and successfully implemented in numerous engineering problems such as the optimum design of truss systems [44,45], steel frames [46], plate girders [47], cylindrical reinforced concrete walls [48], plane stress systems [49], PID controlled active tuned mass dampers [50], retaining walls [51], water network optimization [52], slope stability analysis [53], heat and power systems [54], job shop scheduling [55], team orienteering [56], and vehicle routing [57]. The harmony search method was initially designed for the optimization of discrete-valued variables. Eventually, the method was further developed so that it could be applied to the optimization of continuous-valued data as well. As an example of this kind of optimization problem, the dimensioning of structural components can be given.
To start the harmony search optimization process, the first step is the determination of the design variables and the quantity that is being optimized. In this study, the ply thicknesses (t i ), and the fiber orientation angles (θ i ) of each ply are selected as the design parameters. The quantity that is optimized is the buckling load. The harmony search optimization starts with the random generation of a certain number of design variable combinations. Each one of these design variable combinations is called a candidate solution vector. In the process of generating the candidate solution vectors, certain constraints have to be observed. An example of these constraints in the current study is the minimum ply thickness. Due to constraints related to manufacturing techniques, the layer thicknesses were kept above 0.1 mm in all design vectors. Another design constraint was the total laminate thickness which was not allowed to exceed 2.25 mm to keep the results comparable to the study of Führer [42]. The initial population of candidate solution vectors is updated throughout the harmony search optimization. These updates are in the form of newly generated better-performing solution vectors replacing the vectors generated in the previous steps. The harmony memory consideration rate (HMCR = 0.5 * (1 − iter/maxiter)) and the pitch adjustment rate (PAR = 0.05 * (1 − iter/maxiter)) are the parameters that determine the newly generated solution vectors according to Equation (1). The variables iter and maxiter are the current iteration number and the maximum number of iterations, respectively.
Equation (1) shows the procedure for the calculation of an updated value for the variable with index i in a given candidate solution vector. Here N is the population size, x i,max and x i,min are the maximum and minimum values of the variable with index i in the entire population, respectively, and rand ∈ (0, 1), rand2 ∈ (−1/2, 1/2) are random numbers. In x i,k , i denotes the index of the variable and k denotes the index of the candidate solution vector being used in the update which is chosen as the nearest integer to (rand · N). When the new solution vector generated through Equation (1) performs better than at least one of the existing design vectors, the worst-performing design vector in the population is replaced by the new solution vector. These iterations are repeated until convergence is achieved such that newly generated solution vectors do not perform significantly better than the previous ones. In the case that no further updates are observed after a large number of iterations, the optimization procedure can be aborted. Alternatively, the number of iterations can be used as a stopping criterion for the optimization process. A major limitation of the current study is the limitation on the computing power. Due to limited computing power, the optimization process has to be aborted after a certain number of iterations. Due to the nature of metaheuristic algorithms, the values that the algorithms converge to after a certain number of iterations are not necessarily the exact maximum or minimum values. The harmony search optimization process has been illustrated in Figure 1. The buckling loads of rectangular composite laminates with circular, rectangular, and triangular cut-outs under uniaxial compression are maximized using the harmony search process. Layer thicknesses and fiber orientation angles of the nine layers [42] constitute design vectors with eighteen design parameters. The population of design vectors consists of ten vectors of eighteen variables. The initial population was generated randomly except for the best solution vector which contains the fiber angles and ply thicknesses of the best-performing vector from the previous study of the authors [23] for the corresponding aspect ratio. An illustration of the design variables can be seen in Figure 2. The fiber angle and thickness of a layer are denoted with θ 1 and t 1 , respectively. As it can be seen on the right side of Figure 3, the aspect ratio a/b is the ratio of the longer side length of the rectangular plate to the shorter side length.

Finite Element Model
The plate geometry is discretized with reduced integration elements S4R which is suitable for buckling analysis through linear perturbation that does not demand highperformance computing. Since the linear perturbation analysis computes the bifurcation point of the structure which corresponds to the buckling load, the implementation of geometric imperfections into the structure was omitted. However, to conduct a nonlinear post-buckling analysis it is important to consider the geometric imperfections which can be introduced into the model as a combination of the buckling mode shapes of the structure. In the presence of imperfections, it is expected to observe a smooth transition from the pre-buckling to the post-buckling regime. Furthermore, the addition of reinforcement into the model would give the structure added stiffness and strength which would compensate for the reduced structural stiffness due to the cut-outs. Another favorable effect of the reinforcements is the reduction in the stress concentrations around the cut-out edges. Figure 2 shows one of the simulated plate models with two rectangular cut-outs together with its counterpart without cut-outs to clearly illustrate the differences between the geometries. The used stacking sequence in these models is [ [23]. This combination of fiber angles and plate thicknesses was shown to deliver optimum performance in [23] for rectangular plates of aspect ratio a/b = 2 without cut-outs. In the current study, this stacking sequence has been used as a starting point of the optimization process of plates with cut-outs and it has been further optimized. The boundary conditions are described in Figure 2 as x, z, rz, and MPC where x and z denote the translational constraints in the x and z directions as shown in Figure 2, and rz denotes the rotational constraint around the z-axis. MPC stands for multi-point constraints such that all the nodes on the right edge of the plate are constrained to go through the same amount of displacement as the lower right corner of the plate where a concentrated unit load is applied. This unit load is incrementally increased in the linear perturbation analysis until the critical buckling load is reached [58]. This means that x and z presented the restraint of the translational degrees of freedom at the x and z axes, and rz the restraint of rotational degree of freedom around the z-axis. In addition, all boundary nodes are free to rotate around the x and y axes. To better understand the effect of changing the boundary conditions on the buckling load, a laminated plate with an aspect ratio of two and with circular cut-outs has been simulated with clamped, simply supported, and free support boundary conditions denoted with C, S, and F in Figure 3, respectively. The variation of the buckling load with respect to the boundary conditions was visualized in Figure 3 where the vertical axis shows the normalized buckling load. It can be observed that changing the boundary conditions from clamped to simply supported causes a reduction of 30%, and changing the boundary condition to free support leads to a 50% reduction in the buckling load. Appendix A demonstrates the correlation of the buckling loads obtained using the finite element model with those obtained using analytical techniques.
The material properties of a single layer are listed in Table 1. Here E 1 and E 2 are the Young's moduli of a layer in the fiber direction and perpendicular to the fiber direction, respectively. G 12 and ν 12 are the shear modulus and Poisson's ratio of a lamina, respectively. In ABAQUS, the overall mechanical properties of a laminate are defined by entering the properties of the individual layers constituting the laminate in a sequence. The E 1 and E 2 values combine the properties of the matrix and the fibers in a single value. The harmony search optimization process was controlled by using the PYTHON programming language. A script was programmed to generate the population of solution candidates needed for the harmony search algorithm. For each solution candidate, an ABAQUS input file was generated by using this script, and after the buckling load for the solution candidate was computed with ABAQUS, this buckling load was extracted from the data file generated by ABAQUS using the same PYTHON script. This process is repeated for a predetermined number of harmony search iterations in line with the available computational resources. Table 2 shows a comparison of the buckling loads of the plates with and without circular cut-outs where the optimized stacking sequences from [23] have been used. In Table 2 all buckling loads are obtained through uniaxial compression acting parallel to the longer side of the plate. Clearly, the difference in the buckling load is more pronounced for plates with smaller aspect ratios. In the case of a/b = 1, the reduction in the buckling load due to the cut-out is 19%, whereas in the cases of a/b = 2 and a/b = 3 this reduction is 16% and 8%, respectively. The first buckling modes of a CFRP plate without cut-outs are shown in Figure 4 for all three aspect ratios.  The critical uniaxial compressive load of a symmetric laminated composite plate with simply supported boundary conditions can be calculated using Equation (2), where D 1 , D 2 , D 3 are the plate stiffness values and m, n are the numbers of half-sine waves in the buckled plate in the x and y directions, respectively [59]. Figure 5 conceptually shows the variation of the critical line load for the different aspect ratios when stiffness parameters have fixed values. There are jumps in the value of n x when the mode shape and the number of half-sine waves in x-direction change at a/b = 2 and a/b = 3. After each mode shape changes, the critical load reduces to a level shown with a horizontal dashed line. Figure 5. Variation of the critical distributed load for different aspect ratios.

Results
The first buckling mode shapes for all three aspect ratios are shown in Figure 6 for the circular cut-outs, where the displacement magnitudes are denoted with U. The first buckling mode shapes of plates with rectangular and triangular cut-outs are shown in Figures 7 and 8, respectively. In Figures 4 and 6-8, the buckling mode shapes correspond to uniaxial compression parallel to the longer side of the plate. Figures 6 and 7 compare the difference in the first buckling mode shape caused by changing the cut-out geometry from circular to rectangular. The same difference can also be observed between the mode shapes of the plates with circular and triangular cut-outs.   Table 3 shows a comparison of the buckling loads corresponding to uniaxial compression in the direction of the longer side of the plates for all three aspect ratios and cut-out shapes. Again, the values in this table were obtained using the optimized stacking sequences from Cakiroglu et al. [23]. According to Table 3, the buckling load exhibits a wider variety for the aspect ratio for plates with circular cut-outs. In the case of a/b = 1, the plate with the triangular cut-out delivered the highest buckling load, whereas for a/b = 2 and a/b = 3, plates with circular cut-outs performed better. The variations in the performance of plates with different cut-out geometries can be attributed to the different stress concentrations associated with varying geometries of cut-out. Even though plates with larger a/b values also have greater slenderness, their buckling mode shape has multiple half sinewaves. However, their buckling loads are not necessarily less than plates with smaller aspect ratios. This can also be observed in the values listed in Table 3, where, in the case of circular cut-outs, plates with larger aspect ratios have performed better than their counterparts with smaller aspect ratios.  Figures 9 and 10 are obtained by the harmony search method, the first step of which is the random generation of a population of candidate solutions. The size of this population can be chosen arbitrarily, and in the current study, the population of solution candidates consists of ten members. Each member of the population is represented by a vector containing a total of 18 design variables consisting of 9 layer thicknesses and 9 fiber orientation angles. Once the initial population has been generated randomly, the members of the population are ranked according to their performances. In this process, the best member and the worst member of the population are determined. In the case of buckling analysis, these best and worst members correspond to the members with the highest and lowest buckling load, respectively. After the ranking of the population and the determination of the best and worst members, every member of the population goes through an appropriate harmony search iteration based on Equation (1). This iteration creates new solution vectors, and in case the new solution vectors perform better than the worst member of the population, they replace the worst member. Therefore, after every harmony search iteration, the candidate solution vectors have to be ranked, and the best and worst members of the population have to be determined. Figures 9 and 10 show the performances of these best and worst members of the population throughout the harmony search iterations.   Figure 9 shows the optimization steps for plates without cut-outs for three different aspect ratios. The current study uses the optimized stacking sequences from Figure 9 as a starting point for the optimization of plates with cut-outs. Figures 10-12 show the optimization steps for plates with circular, rectangular, and triangular cut-outs, respectively where the black curves show the development of the best solution, and the red curves show the development of the worst solution in the population of solution candidates. The outcome of this optimization is listed in Tables 4-6 for circular, rectangular, and triangular cut-outs, respectively. It should be noted that all the values listed in Tables 4-6 are the result of uniaxial compression parallel to the longer side of the plate. Table 4 shows that, for all aspect ratios, the optimized stacking sequences performed better than the conventional stacking sequence of [42] in the case of circular cut-outs. As Table 5 shows, the optimized stacking sequences performed better than the conventional symmetric sequences in the case of rectangular cut-outs for the aspect ratios of a/b = 1 and a/b = 3. Particularly in the case of a/b = 3, the increase in the buckling load was 4.6%. For the other two aspect ratios, the differences between the optimized and symmetric sequences were less than 2.7% in the case of rectangular cut-outs. Similarly, for the case of triangular cut-outs, the optimized stacking sequences and their performance comparisons can be found in Table 6. In this case, the optimized stacking sequence performed 1.5% better than the symmetric sequence for the aspect ratio of a/b = 2. For the other two aspect ratios, the symmetric stacking sequence performed better than the optimized sequence. A summary of the results from Tables 4-6 is presented in Figure 13, where the variation of the buckling load (P cr ) for the cut-out shape, and the aspect ratio is plotted. Due to the availability of computational resources, all optimizations were aborted after 1200 harmony search iterations which is a limitation of the current study.     Figure 13. Variation of the buckling load with respect to cut-out shape and aspect ratio.

Discussion
Carbon fiber reinforced laminated composite plates of three different aspect ratios are simulated in the presence of circular, square, and triangular cut-outs. For all three cut-out shapes, the stacking sequence has been optimized using the harmony search algorithm. As the initial configuration that is being improved, the best stacking sequence obtained by the authors [23] has been used. Therefore, as shown in Figures 10-12, the improvement steps in the best solution vector are sparse and of relatively small magnitude. On the other hand, the improvement steps of the worst solution vectors are significantly larger in the beginning. They get smaller eventually because these vectors have been initialized randomly, rather than the best solution vector, which was initialized with the values of an already optimized solution vector.
Using the optimized stacking sequence from Cakiroglu et al. [23], the buckling loads are computed and listed in Table 3 for all aspect ratios and cut-out geometries. Table 3 shows that the increase in the aspect ratio did not have a considerable effect on the buckling load for square and triangular cut-outs while it increased the buckling load of the plates with circular cut-outs. On the other hand, the buckling load values obtained through the symmetric stacking sequence of [42] and listed in Tables 4-6 show that an increase in the aspect ratio is accompanied by a slight increase in the buckling load for all cut-out geometries. It can also be observed that, according to Tables 4-6, for all aspect ratios, the highest value of the buckling load was obtained from the plate with the triangular cut-out.
The optimization results in Table 4 clearly show that laminates with cut-outs with optimized stacking sequences perform better than the regular stacking sequence identified by Führer [42] for all aspect ratios. According to Table 4, the optimized laminates performed up to 3.43% better than their regular counterparts. The best performance improvement was observed for a/b = 3 with a fiber angle sequence of [ where the angle is in degrees and the ply thickness is in mm. According to Table 5, the optimization algorithm obtained a stacking sequence that performs better than the symmetric stacking sequence for the aspect ratios a/b = 1 and a/b = 3 in the case of a rectangular cut-out. The best performance improvement was observed for the angle sequence of [ Further research in this area can be carried out on the effects of changing the size and location of the cut-outs. The effects of using fibers with curvilinear paths on the performance of the laminate could also be investigated. Another direction for future research could be the optimization of laminates with cut-outs made of different composite materials such as boron/epoxy or glass/epoxy composites.

Conclusions
Cut-outs in laminated composite plates are prevalent, and they play a crucial role in the performance of a structure. Since these structures are primarily thin-walled, they are prone to buckling, and maximizing their buckling loads bears great significance. Previously obtained solutions have been used as a starting point in the current study, and they have been further optimized with the addition of circular, rectangular and triangular cut-outs into the finite element models. The effect of changing the cut-out geometry on the buckling load has been demonstrated. As the optimization methodology, the well-established harmony search algorithm has been applied. The key outcomes of this study can be summarized as follows:

•
The application of dispersed stacking sequences where the fiber orientation angles can take any value between 0 • and 90 • , can increase the buckling load of laminated composite plates. For three different aspect ratios, the optimized stacking sequences have been found to perform better than their counterparts with regular stacking sequences consisting of 0 • , 45 • , and 90 • fiber orientation angles. • On average, plates with triangular cut-outs performed better than plates with circular and rectangular cut-outs, while only minor performance differences could be observed between plates with rectangular and circular cut-outs.

•
The harmony search algorithm could be effectively applied to the problem of buckling load maximization for plates with cut-outs. Improvements in the structural performance could be achieved with relatively minor computational effort. • For all aspect ratios and cut-out geometries, a lack of a regular pattern could be observed in the optimized sequences of fiber angles and layer thicknesses which necessitates the application of optimization algorithms such as the harmony search algorithm on a case-by-case basis.
Author Contributions: C.C. and G.B. generated the analysis codes; the text of the paper was formed by C.C., K.I., G.B., S.K. and Z.W.G.; the figures were prepared by C.C. and K.I.; supervised the research direction, G.B. and Z.W.G. All authors have read and agreed to the published version of the manuscript.

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

Verification of the Finite Element Model by Analytical Techniques
In the case of a square plate with a symmetric stacking sequence, the equation for the buckling load of a plate is given in Equation (A1) [59]. In Equation (A1) n cr is the critical free edge load acting on one side of the plate, b is the side length of the square plate and D 1 , D 2 , D 3 are the components of the bending stiffness matrix D ij with D 1 = D 11 , D 2 = D 22 , and D 3 = D 12 + 2D 33 . The bending stiffness matrix is related to the lamina stiffness matrices Q and Q given in Equation (A2) where T is the transformation matrix and θ is the orientation angle of a lamina. The components of the bending stiffness matrix can be calculated as in Equation (A3) where z k defines the distance of a lamina from the mid-plane of the laminate as shown in Figure A1 and N is the total number of layers. Furthermore, the fiber orientation angle θ and its position with respect to the global xy-coordinates is shown in Figure A2.  The total thickness and side length of the plate are chosen to be 2.25 mm and 150 mm, respectively, where all layers have the same thickness. A total of 12 stacking sequences with fiber orientation angles of 0 • , 45 • , 60 • and 90 • have been analyzed. Figure A3 shows the correlation between the analytical results and the finite element solutions. A coefficient of determination of R 2 = 0.96 could be observed.