A novel multi-objective optimization framework for tensile properties of carbon/Kevlar fiber-reinforced composite laminates based on improved hybrid SPEA2-MOPSO algorithm and TOPSIS method

Abstract

Hybrid fiber-reinforced composite laminates are extensively employed across various engineering applications owing to their superior comprehensive mechanical performance. However, effectively balancing competing performance parameters remains critical to maximizing their overall structural capabilities and meeting practical engineering requirements. To address this challenge, this study presents a multi-objective optimization framework for the design of carbon/Kevlar fiber-reinforced composite laminates, with a focus on the influence of ply material and stacking sequence on tensile behavior under load. Initially, optimal Latin hypercube sampling (LHS) was employed to generate training samples, and a corresponding simulation database was established via finite element analysis to construct a surrogate model. A novel optimization framework was then developed by integrating a Kriging surrogate model with a hybrid strength Pareto evolutionary algorithm 2–multi-objective particle swarm optimization (SPEA2-MOPSO) algorithm and the technique for order preference by similarity to ideal solution (TOPSIS). Within this framework, the ply sequence and material type served as decision variables, while tensile strength, tensile modulus, and elongation at break were considered as optimization objectives. The research results show that the finite element simulation using the progressive damage model was verified by the tensile test, which confirmed the accuracy of damage prediction and revealed the synergistic effect among tensile performance parameters. A comparative analysis of three optimization algorithms, namely MOPSO, SPEA2, and the hybrid SPEA2-MOPSO, demonstrated the superior accuracy and efficiency of the hybrid approach in solving complex multi-objective problems. Finally, the application of the TOPSIS method enabled the identification of optimal trade-off solutions tailored to various engineering priorities, providing practical guidance for structural design decisions.

Keywords

hybrid fiber composites kriging surrogate model hybrid SPEA2-MOPSO algorithm multi-objective optimization

Introduction

Carbon fiber-reinforced polymers (CFRPs) are widely utilized as lightweight structural materials in industrial applications due to their high specific strength and stiffness. However, their widespread adoption is often limited by high production costs, intrinsic brittleness, and susceptibility to catastrophic fracture. In comparison, Kevlar fiber-reinforced polymers (KFRPs) offer improved toughness and cost-effectiveness, making them attractive alternatives. To enhance the mechanical performance of hybrid fiber-reinforced polymer (HFRP) composite laminates and to achieve an optimal balance between weight reduction and economic feasibility, relevant research has received growing attention in recent years. Wang et al. (2024) studied the mechanical properties of carbon/ Kevlar fiber composite laminates by combining carbon fibers and Kevlar fibers in the same matrix. The results demonstrated that the hybrid composites can effectively leverage the unique advantages of each constituent material. With the continuous advancement of research and development in composite materials, their properties and performance have been improved, and their future application prospects will become more extensive, particularly in wind energy, automotive, rail, light aviation, construction, and other key industries (Berton et al., 2026; Chairi et al., 2026; Hoshikawa et al., 2025; Slamani and Chatelain, 2015; Wang et al., 2023). Although HFRP composite laminates demonstrate favorable tensile properties, a comprehensive understanding of their overall tensile behavior is essential for optimizing performance and ensuring their reliability in demanding engineering applications.

HFRP composite laminates are multilayered structures composed of various fiber types embedded within a polymer matrix, arranged in specific stacking sequences and layer configurations (Peng et al., 2025). As a result, they demonstrate distinctive mechanical characteristics, including high tensile strength, specific strength, specific modulus, elongation at fracture, long fatigue life, and excellent corrosion resistance. However, significant variations in their performance can arise due to differences in the constituent materials and the adopted lay-up configurations. Zheng et al. (2017) and Wang et al. (2022) both improved the mechanical properties of composites by incorporating aramid fibers (Kevlar-fiber) with higher tenacity. Their findings indicate that differences in lay-up configurations significantly affect the tensile strength and toughness of HFRP composites. With the advancement of computational technologies in recent years, finite element method (FEM)-based numerical simulation has emerged as a cost-effective and widely used approach for analyzing the mechanical behavior of composite materials (Houari et al., 2026; Zhang et al., 2025). Under tensile loading, composite laminates commonly exhibit dimensional fracture, matrix cracking, and interfacial debonding. In order to most accurately simulate the damage behavior of laminates, based on specific damage theories and assumptions, scholars have proposed many numerical modeling methods using different failure criteria and damage evolution methods (Zhou et al., 2020). Due to the differences in numerical modeling methods under different working conditions, adopting appropriate failure criteria and evolution methods can effectively improve the accuracy of numerical simulations. Li et al. (2019) studied composite laminates based on different failure criteria and damage evolution methods to compare their applicability. Their results showed that while the Hashin criterion was inadequate for simulating matrix tensile failure, it was more effective in capturing compressive damage. Similarly, Huang et al. (2023) investigated failure criteria, evolution models, and interface modeling techniques. Their results demonstrated that the combination of the Puck criterion, linear damage evolution, and cohesive zone elements (with either finite or zero thickness) produced simulations more consistent with experimental data. Through the FEM, these studies aim to determine the configuration with the best mechanical properties among a selected finite number of laminate structures. While FEM provides valuable insights and enables the selection of laminate configurations with superior mechanical performance from a predefined set of options, its reliance on a finite sample space limits the potential to identify globally optimal solutions (Zhang et al., 2025).

To optimize the performance of composite laminates, genetic algorithms (GAs) have gained significant attention due to their robustness and ability to efficiently explore global solution spaces. Compared to finite element simulations, which are often computationally expensive and limited by the number of allowable simulation runs, GAs offer a more scalable approach to structural optimization. Consequently, in recent years, GAs have been widely adopted for their effectiveness in identifying global optima in complex structural systems (Chen et al., 2023). However, directly coupling GAs with FE simulations substantially increases computational time and cost, thereby reducing overall optimization efficiency. To address this issue, researchers have introduced surrogate modeling techniques, such as response surface methodology (RSM) (Hou et al., 2023). These surrogate models are constructed using a limited number of experimental and simulation data points, enabling efficient function approximation. When integrated with GAs, surrogate models facilitate fast and accurate optimization, striking a balance between computational cost and solution accuracy. As a result, the hybridization of GAs and surrogate models has become increasingly popular in recent composite optimization research (Ghane et al., 2024; Sahu and Andhare, 2019). For example, Mostofi et al. (2020) and (Arhore and Yasaee, 2020) carried out the structural optimization design of laminates by combining the response surface method with the GA, and verified the effectiveness of the optimization scheme. These studies validate the potential of coupling surrogate models with evolutionary algorithms to address the trade-off between simulation cost and performance precision. Furthermore, the optimization of hybrid fiber-reinforced composite laminates constitutes a complex multi-objective optimization problem involving several discrete variables. Key design parameters typically include stacking material (Gholami et al., 2021), stacking sequence (Li et al., 2019), layer thickness (Garg et al., 2023), and fiber orientation angle (Garg et al., 2023), among others. Numerous researchers have explored the optimization of these variables. For instance, (Beylergil, 2020) used a multi-objective genetic algorithm (MOGA) to conduct design optimization for the weight and stiffness of hybrid composite laminates subjected to eccentric loads, and determined the optimal lay-up sequence and lay-up materials (carbon fibers or E-glass). Similarly, Fakoor et al. (2019) optimized the weight and buckling load of hybrid laminated composite plates through the MOGA, regarding the orientation, thickness, quantity, and material of each layer as design variables. Although these approaches have successfully produced optimized solutions, they are generally well-suited for problems with low-dimensional design spaces. In high-dimensional scenarios, the efficiency of GAs tends to decline, increasing the likelihood of premature convergence and the inability to locate the global optimum (Ghiasi et al., 2009). Therefore, it is essential to develop effective dimensionality reduction techniques for representing lay-up configurations and to design optimization algorithms with improved exploration capabilities to overcome convergence limitations in high-dimensional composite design problems.

Based on previous studies, it has been observed that some researchers have explored dimensionality reduction techniques for the lay-up parameters of single-fiber composite laminates (Imran et al., 2021; Xu and Wang, 2022). However, the design and optimization of HFRP composite laminates present greater complexity, and related research remains relatively limited. For instance, from the perspective of optimization algorithms, there is a need to balance optimization accuracy with convergence speed, which is a trade-off that is difficult to achieve using traditional single algorithms alone (Xu et al., 2024). As a result, hybrid algorithmic frameworks, which integrate multiple optimization strategies, have been developed to enhance performance in such scenarios. Wang et al. (Wang and Cai, 2018) proposed a hybrid method that couples the particle swarm optimization (PSO) algorithm with the bacterial foraging optimization (BFO) algorithm to improve crashworthiness and lightweight design. This hybrid algorithm reduces computational time in obtaining the Pareto solution set and avoids the optimization from falling into local optimality. Mir and Rezaeian (2016) proposed an effective method based on a hybrid of the PSO and the GA to solve the parallel machine scheduling problem. This method proved both effective and feasible in generating optimal or near-optimal solutions, especially under conditions of increased problem complexity. Comparably, Tsai et al. (2004) proposed a hybrid Taguchi's method and GA (HTGA) to solve multi-objective optimization problems with continuous variables. The results show that not only can optimal or near-optimal solutions be found, but also better robustness results can be obtained. Despite these advancements, the integration of hybrid optimization algorithms specifically tailored for the design optimization of HFRP laminates has not been thoroughly investigated. Further research is needed to enhance algorithm applicability and efficiency in this domain.

Building upon previous research and existing methodologies, this study focuses on the mechanical performance and design flexibility of HFRP laminates. A progressive tensile damage simulation model is developed for carbon/Kevlar fiber hybrid composite laminates to investigate the effects of stacking sequence and hybridization ratio (i.e., the proportion of different fiber materials) on their tensile behavior and failure mechanisms. To facilitate structural representation, a real-number encoding scheme is proposed to describe various combinations of lay-up materials and stacking sequences. A multi-objective optimization framework is constructed by integrating a Kriging surrogate model with the strength Pareto evolutionary algorithm 2 (SPEA2) and multi-objective particle swarm optimization (MOPSO), aimed at enhancing the overall tensile performance of the hybrid laminates. Furthermore, to support optimal solution selection from the obtained Pareto front, the technique for order preference by similarity to ideal solution (TOPSIS) is employed in conjunction with principal component analysis (PCA) to assign weights to multiple objectives. This approach allows for flexible decision-making based on varying design preferences and performance requirements.

Materials and methods

Material details

The raw materials used in the fabrication of HFRP composite laminates include plain weave carbon fiber fabric, plain weave aramid (Kevlar) fiber fabric, epoxy resin, and a curing agent. The material properties and corresponding manufacturers are detailed in Table 1.

Table 1.

Material characteristics and manufacturer of raw materials.

Raw material	Material characteristic	Manufacturer
Carbon fiber fabric	Thickness:0.111 mm, areal density: 200 g/m², fiber type: 3 K, plain-weave	Zhongfu Shenying Carbon fiber Co., Ltd., China
Aramid fiber fabric	Thickness:0.127mm areal density: 200 g/m², fibertype:Kevlar-29,1500D, plain-weave	Dingyuan Carbon Fiber Technology Co., Ltd., China
Epoxy resin	105 epoxy resin, toxic, boiling point>204 $^{\circ} C$	West System, Inc., USA
Hardener	206 slow hardener, toxic, boiling point >204 $^{\circ} C$	West System, Inc., USA

Material failure theory

Damage models of intra- and inter-laminar

The primary failure mechanisms in composite materials include fiber tensile failure, fiber compressive failure, matrix cracking, and matrix fracture (Wang et al., 2024). Under tensile loading, failure in HFRP laminates typically results from a combination of fiber tensile failure, matrix cracking, and matrix fracture. This indicates the presence of a multi-mode failure mechanism. Commonly used failure criteria for modeling such behavior include the Hashin, Chang–Chang, and Hou failure criteria. To accurately simulate the progressive failure behavior of HFRP laminates under tensile loading, a progressive damage model based on the Chang–Chang failure criterion (Material Model 54/55—Enhanced Composite Damage Model) was implemented using HyperMesh and LS-DYNA. This model incorporates a bilinear constitutive law to simulate the initiation and propagation of interfacial damage under mixed-mode conditions. Based on previous research findings (Jiang et al., 2020), the intra-laminar and inter-laminar damage parameters used in this study are summarized in Table 2.

Table 2.

The intra- and inter-laminar damage models.

Intralaminar model
Failure modes	Failure criterion	Strength reduction
Tensile fiber failure	$e_{f}^{2} = {(\frac{σ_{a a}}{X_{t}})}^{2} + β {(\frac{σ_{a b}}{S_{c}})}^{2}$	$X_{t} = 0.07 X_{t}$ (1)
Compressive fiber failure	$e_{c}^{2} = {(\frac{σ_{a a}}{X_{c}})}^{2}$	$X_{c} = 0.14 X_{c}$ (2)
Tensile matrix failure	$e_{m}^{2} = {(\frac{σ_{b b}}{Y_{t}})}^{2} + {(\frac{σ_{a b}}{S_{c}})}^{2}$	$Y_{t} = 0.2 Y_{t}$ (3)
Compressive matrix failure	$e_{d}^{2} = {(\frac{σ_{b b}}{2 S_{c}})}^{2} + [{(\frac{Y_{c}}{2 S_{c}} - 1)}^{2}] \frac{σ_{b b}}{Y_{c}} + {(\frac{σ_{a b}}{S_{c}})}^{2}$	$X_{c} = 0.4 X_{c}$ (4)
Interlaminar and interface model
Failure criterion	Damage evolution criterion	Mix-mode fracture energy criterion
${(\frac{t_{n}}{t_{n}^{0}})}^{2} + {(\frac{t_{s}}{t_{s}^{0}})}^{2} + {(\frac{t_{t}}{t_{t}^{0}})}^{2} \leq 1$ (5)	$D = \frac{δ_{m}^{f} (δ_{m}^{m a x} - δ_{m}^{0})}{δ_{m}^{m a x} (δ_{m}^{f} - δ_{m}^{0})}$ (6)	$G^{C} = G_{n}^{C} + (G_{s}^{C} - G_{n}^{C}) {(\frac{G_{s}}{G_{T}})}^{η}$ (7)

Among them, the formulas for the four failure modes of the intra-layer model are Eqs. (1)–(4). Where $e_{f}$ , $e_{c}$ , $e_{m}$ , $e_{d}$ are damage factors corresponding to four failure modes. $σ_{a a}$ , $σ_{a b}$ , $σ_{b b}$ are the effective stress. $X_{t}$ , $X_{c}$ , $Y_{t}$ , $Y_{c}$ , $S_{c}$ are fiber, matrix, and shear strengths, respectively.

The interlaminar damage of the HFRP laminate is simulated using cohesive elements, and the damage initiation point of the cohesive elements is judged by the quadratic stress criterion in Eq. (5). Where $t_{n}$ , $t_{s}$ , and $t_{t}$ represent the normal and shear tractions; $t_{n}^{0}$ , $t_{s}^{0}$ , and $t_{t}^{0}$ denote the interface normal and shear strengths.

The damage evolution criterion of cohesion adopts the energy-based mixed-mode B-K criterion, and its expression is Eq. (7). Where $G^{C}$ is the fracture toughness in mixed mode. $G_{n}^{C}$ , $G_{s}^{C}$ are the critical strain energy release rates in tangential and normal directions. $G_{S}$ , $G_{T}$ are the energy release rates corresponding to the damage modes of the cohesive element and $η$ is the material parameter.

Table 3 presents the numerical simulation parameters for the properties of CFRP, KFRP, and interlaminar cohesive elements. Where $G_{f}^{1 T}$ , $G_{f}^{1 C}$ , $G_{f}^{2 T}$ , and $G_{f}^{2 C}$ are the longitudinal tensile and compressive, transverse tensile and compressive fracture toughness respectively. $E_{N}$ and $E_{T}$ are respectively the stiffness in the plane of the cohesive element and normal to its plane. $T$ and S are respectively the peak traction in the normal direction and the peak traction in the tangential direction. $G_{IC}$ and $G_{IIC}$ are the energy release rates for mode I and mode II, respectively.

Table 3.

Material parameters used for numerical simulations, including the intralaminar properties of HFRP (Bresciani et al., 2016; Gower et al., 2008; Sokolinsky, 2011; Zhao et al., 2019) and the parameters of the interlaminar cohesive element (Xu et al., 2024).

	Intralaminar properties of HFRP
parameter	CFRP	KFRP	parameter	CFRP	KFRP
$ρ /$ g·cm⁻³	1.44	1.23	$X_{T}$ /MPa	579.67	255
$E_{11}$ /GPa	57.02	18.5	$X_{C}$ /MPa	412.53	185
$E_{22}$ = $E_{33}$ /GPa	57.02	18.5	$Y_{T}$ /MPa	579.67	255
$G_{12}$ = $G_{13}$ /GPa	4.586	7.7	$Y_{C}$ /MPa	412.53	185
$G_{23}$ /GPa	4.586	7.7	$S$ /MPa	68	56.7
$v_{12}$ = $v_{13}$	0.067	0.25	$v_{23}$	0.032	0.33
$G_{f}^{1 T}$ /kJ·m⁻²	85.4	170.8	$G_{f}^{2 T}$ /kJ·m⁻²	85.4	170.8
$G_{f}^{1 C}$ /kJ·m⁻²	155.3	310.6	$G_{f}^{2 C}$ /kJ·m⁻²	155.3	310.6
	Interlaminar cohesive element
$ρ /$ g·cm⁻³	E_N/GPa	E_T/GPa	T/MPa	S/MPa	G_IC =G_IIC/J·m⁻²
1.2	30	30	36	36	0.25

Methods

Hybrid structure design

To investigate the effect of hybrid fiber composite laminate architecture on tensile properties, this study focused on carbon and aramid fibers by designing seven composite laminate configurations, each comprising nine layers [Fig. 1(b)]. These configurations included different hybridization ratios, namely CF8/KF1, CF5/KF4, and CF2/KF7, as well as varied stacking sequences, specifically CF5/KF4, KF2/CF5/KF2, and CF/(KF/CF)_2s. Additionally, two monolithic configurations, pure carbon fiber (CF9) and pure aramid fiber (KF9) laminates, were fabricated for comparative analysis.

Figure 1.

Tensile specimen design and testing. (a) Test machine. (b) Schematic diagrams of HFRP design. (c) Tensile specimen size.

Fabrication of samples

Based on the proposed structural designs, all laminates were manufactured using woven carbon and Kevlar fiber fabrics as reinforcements. The matrix material consisted of 105 epoxy resin and 206 slow hardener, mixed at an optimal ratio of 5:1 by weight at room temperature (Cheng et al., 2020). Laminate fabrication was performed using a hydraulic press at 1.5 MPa for 12 h. Following depressurization, the semi-cured laminates were post-cured in an air oven at 60 °C for 3–4 days (Jiang et al., 2021; Jiang and Ren, 2022). Final laminate dimensions were 110 mm × 20 mm × 1.8 mm.

Experiment and numerical simulation

Tensile properties test

Tensile tests were conducted according to the ASTM D3039 standard (ASTM D3039, n.d.) using a CTM8010 microcomputer-controlled electronic universal testing machine equipped with a 20 kN load cell [Fig 1(a)]. One end of each specimen was clamped, while the other end was pulled at a constant displacement rate of 2 mm/min. The test was terminated automatically upon fracture, indicated by a sudden drop in tensile force. To prevent premature failure due to stress concentrations at the clamping ends, aluminum tabs (25 mm × 20 mm × 2 mm) were bonded to both ends of the specimen using epoxy adhesive [Fig. 1(c)]. These tabs provided a uniform gripping surface, transferred load efficiently, and protected the outer fibers from damage (Wang et al., 2025).

To ensure reliability of the experimental results, five specimens from each test case were loaded until failure and the average value was used as the final result. After testing, the stress–strain data were extracted, and the tensile modulus and elongation at break were calculated using Eqs. (8) and (9). The resulting tensile properties are summarized in Table 4.

E_{a} = \frac{Δ σ_{t}}{Δ ε_{L}}

(8)

ε = \frac{Δ L}{L} \times 100 %

(9)

Table 4.

Tensile properties of composite laminates.

Hybrid sample	Tensile strength	Tensile modulus	elongation at peak
Hybrid sample	(MPa)	(GPa)	(%)
CF9	628.65	232.22	5.81
CF/KF-HFRP
CF8/KF1	619.13	183.44	5.85
CF5/KF4	536.49	178.00	6.41
CF2/KF7	466.86	172.70	7.16
KF9	521.23	135.21	8.42
Stacking sequence
CF5/KF4	536.49	178.00	6.41
KF2/CF5/KF2	482.66	160.42	5.32
CF/(KF/CF)2s	440.30	156.65	4.78

Where $E_{a}$ is the apparent tensile modulus, $Δ σ_{t} / Δ ε_{L}$ is the slope of the initial linear portion of the plotted stress–strain curve. $ε$ is the elongation at break， L the length of the sample， $Δ L$ the elongation corresponding to the break of the sample.

Finite element simulation and verification

Finite element model (FEM) and simulation

The finite element model of the tensile test is shown in Fig. 2(a). HFRP laminates were modeled using the MAT_054 material model in LS-DYNA, implemented via HyperMesh. A progressive failure model based on the modified Chang–Chang failure criterion was used to simulate fiber and matrix failure, with material properties defined as per Table 2. The laminate was constructed with nine layers, each 0.2-mm thick. A mesh convergence study [Fig. 2(c)] was performed to ensure mesh quality. The model employed 1 mm × 1 mm hexahedral solid elements for in-plane meshing. To simulate interlaminar failure, zero-thickness eight-noded cohesive elements were inserted between each pair of layers [Fig. 2(b)]. AUTOMATIC_SINGLE_SURFACE contact was established between the cohesive elements and adjacent layers, with a tangential and normal friction coefficient of 0.2 employed. Each ply was assigned its respective material properties and lay-up angle through a local coordinate system: 0° aligned with the global X-axis, 90° with the Y-axis, and the Z-axis representing the through-thickness direction.

Figure 2.

Tensile simulation model of the hybrid laminate. (a) Finite element model. (b) Stacking structure and cohesive elements. (c) Mesh independence analysis.

Boundary conditions were applied to simulate clamped ends during tension. Rigid constraints were imposed at both ends of the laminate [Fig. 2(a)] to mimic the non-movable, fixed boundary during testing. All degrees of freedom (DOFs) were constrained at one end of the specimen (Azhdari and Taheri-Behrooz, 2023), while the opposite end was allowed to move only along the X-axis at a constant speed corresponding to the experimental tensile rate.

Model verification

The simulation results obtained using HyperMesh-LSDYNA and LSPP software were compared with experimental data, as illustrated in Fig. 3. The stress–strain responses from the simulations demonstrate good agreement with the experimental curves, thereby validating the accuracy and feasibility of the developed finite element model. Fig. 3(a) presents the stress–strain curves of laminates with different hybridization ratios under both experimental and simulation conditions. As the proportion of carbon fiber layers decreases (i.e., the proportion of aramid fiber increases), a gradual reduction in the maximum tensile strength is observed. Similarly, Fig. 3(b) shows the tensile response of laminates with different stacking sequences. It is evident that as the number of carbon/Kevlar fiber (CF/KF) interfaces increases, the maximum tensile strength decreases accordingly. Fig. 3(c) and 3(d) illustrates the comparison of tensile modulus and elongation at break, respectively, across various laminate configurations under both experimental and simulated conditions. The tensile modulus decreases with an increasing KF layers and with the number of CF/KF interfaces. The relative error between experimental and simulated tensile modulus ranges from a minimum of 0.31% to a maximum of 6.48%. Conversely, the elongation at break increases with a higher hybridization ratio (i.e., more KF layers) but decreases as the number of CF/KF interfaces rises. The corresponding relative error ranges from 1.92% to 7.18%. These results further confirm that both the stacking sequence and hybridization ratio have a significant impact on the tensile behavior of hybrid laminates. Finally, it can also be concluded that among all these designed laminates, CF5/KF4 has more outstanding comprehensive tensile properties.

Figure 3.

Comparison of tensile property results between experimental tests and finite element model. (a) The different hybrid ratios. (b) The different stacking sequences. (c) The tensile modulus. (d) The elongation at peak.

Fig. 4 compares the failure morphology of laminates with different hybridization ratios and stacking sequences under tensile loading. It is noteworthy that tensile failure is accompanied by fiber breakage, fiber pull-out, and interlaminar delamination—phenomena consistent with the findings of Wang et al. (2025). Their work identified common tensile failure modes in composite laminates, including fiber breakage, pull-out, fiber/matrix debonding, matrix fracture, delamination, and interfacial failure. In laminates with fewer KF layers [Fig. 4(a)], the fracture path is visible, but delamination is not pronounced. As the KF layer increases, the laminates exhibit significant ductile failure behavior characteristic of Kevlar fibers. This increased ductility induces fiber buckling and causes carbon fiber layers to warp upward, thereby intensifying delamination [Fig. 4(b)]. Fig. 4(c) further shows that even when the carbon fiber layers have fractured, the laminate retains some load-bearing capacity due to the continued contribution of pulled-out Kevlar fiber bundles, delaying complete failure. These observed failure mechanisms demonstrate that variations in hybridization ratio introduce multiple failure modes, leading to increasingly complex fracture patterns under tensile loading. In particular, delamination predominantly occurs at CF/KF interfaces. As shown in Fig. 4(e), laminates with a higher number of layer interfaces experience more frequent and severe interfacial delaminations. Therefore, increasing interface quantity due to altered stacking sequences exacerbates damage complexity and contributes to a reduction in overall laminate performance.

Figure 4.

Comparison between experiment and numerical simulation on tensile failure behavior and damage appearance. (a) The CF8/KF1 laminate. (b) The CF5/KF4 laminate. (c) The CF2/KF7 laminate. (d) The KF2/CF5/KF2 laminate. (e) The CF/(KF/CF)2s laminate.

Laminated composite optimization and multi-objective algorithm

Coding strategies for laminated composite optimization

Based on the experimental and simulation analyses presented in the preceding sections, it can be concluded that both ply material and stacking sequence have a significant influence on the tensile properties of laminated composite plates. Due to the synergistic interaction between ply material and stacking sequence (Wang et al., 2025), independently optimizing these two parameters may fail to achieve a globally optimal solution. To address this, the present study proposes an integrated optimization strategy in which the ply material and stacking sequence of each individual layer are treated as discrete design variables and optimized simultaneously. The candidate ply materials are carbon fiber and Kevlar fiber, with the total number of laminate layers fixed at nine. Given that the design variables are discrete in nature, a real-number encoding scheme is adopted to facilitate multi-objective optimization. Specifically, the real numbers ‘1’ and ‘2’ are used to represent CF and KF, respectively. An illustration of this encoding strategy is provided in Fig. 5, which visually demonstrates how the stacking sequence and material assignment are digitally implemented. This real-number encoding framework is programed using MATLAB, allowing for the rapid generation of laminate configurations during the optimization process. By encoding both ply type and sequence into a unified design representation, the proposed method enables a more efficient and comprehensive search of the solution space.

Figure 5.

Stacking sequence coding strategy.

Definition of the optimization problem

As analyzed in Section 3, a positive correlation exists between the tensile strength and tensile modulus of laminated plates, which means that higher tensile strength is generally accompanied by a higher tensile modulus. However, this increase in modulus also leads to elevated overall stiffness, which in turn reduces the elongation at break. These observations indicate a coupling effect among the three mechanical properties. To address this issue, a multi-objective design approach is required to identify an optimal trade-off solution that balances these competing factors. In this study, tensile strength and elongation at break are selected as the primary optimization objectives, as they represent the most critical performance indicators for structural integrity and flexibility. Meanwhile, tensile modulus and a cost index are incorporated as constraint conditions to ensure practical feasibility and economic efficiency. The mathematical formulation of the multi-objective optimization model is expressed as follows:

f i n d X = [x_{1}, x_{2}, x_{3}, \cdot \cdot \cdot, x_{9}]

x_{i} \in 1, 2

max f (x) = f_{S} (x), f_{M} (x), f_{E} (x)

(10)

s . t . {\begin{matrix} f_{M} (x) > 135.21 \\ f_{C} (x) < 0.37 \end{matrix}

where x represents the design variables for the lay-up materials and lay-up methods. $f_{S} (x)$ and $f_{E} (x)$ represent the functions of tensile strength and elongation at break, respectively, while $f_{M} (x)$ and $f_{C} (x)$ represent the functions of tensile modulus and laminated plate cost, respectively. The cost of carbon fiber raw material is 21.07 $/m², and the cost of Kevlar-29 fiber fabric raw material is 15.22 $/m². The area of a single layer is 0.0022 m².

Kriging surrogate model

The Latin hypercube sampling (LHS) method (Zhang et al., 2022) was employed to generate 100 sets of sample points, and the corresponding tensile properties were obtained through finite element numerical simulations. Table 5 presents a subset of the simulation results, while Fig. 6 illustrates the data distribution using violin plots, providing insight into the variability and density of the tensile performance indicators. Subsequently, the Kriging surrogate model (Lu et al., 2023), implemented in MATLAB, was utilized to establish the functional relationship between the design variables and the tensile performance responses.

Figure 6.

Violin plots for the numerical simulation data. (a) The tensile strength. (b) The tensile modulus. (c) The elongation at peak.

Table 5.

Tensile response values of composite laminates. (C represents the carbon fiber layer, and K represents the Kevlar fiber layer).

FEMnumber	Stackingstructure	Tensile strength(MPa)	Tensile modulus(GPa)	Elongation at peak(%)
1	CCCCKCCCC	571.70	182.42	5.78
2	CCCCCCKCC	571.37	182.21	5.72
3	CCCCCCCKC	581.16	181.37	6.04
4	CCCCCCCKK	617.94	182.94	5.82
5	KCCCCCCCK	546.47	175.36	5.81
6	CCKCCKCCC	513.00	171.04	5.51
7	CCKCCCKCC	536.49	174.24	5.67
8	CCCCKKCCC	533.56	173.68	5.63
9	KKCCCKKCC	468.62	158.28	5.22
10	CCKKCKKCC	456.39	157.93	5.15
11	CCCCCCKKK	565.30	179.30	6.22
12	CCCKKKCCC	531.00	172.56	6.04
13	CCCCKKCCK	521.25	169.18	5.99
14	CCKKCCKCC	440.98	162.44	5.91
15	CCKKKKKCC	519.85	174.61	6.36
16	KKCCKCCKK	501.93	157.59	6.13
17	CKCKKKCKC	433.87	156.47	5.39
18	KCKCKCKCK	400.56	155.35	5.01
19	CCCKKKKKK	478.35	173.10	6.78
20	KKKCCCKKK	477.04	155.52	6.55
21	CKKKKKKKC	457.73	154.56	6.79
22	KKKCCKKKK	458.76	155.12	6.93
23	KCKKKCKKK	455.53	144.56	6.47
24	KCKKKKKCK	421.25	143.44	6.19
25	KKKCKCKKK	451.99	144.00	6.33
26	CKKKKKKKK	459.89	170.63	7.54
27	KKKKCKKKK	443.47	151.37	7.17
28	KKKKKCKKK	449.94	151.93	7.31
29	KKKKKKCKK	434.42	152.49	7.03
30	KKKKKKKCK	452.27	153.05	6.89

An additional 10 sets of validation sample points and five sets of experimental test values are provided in Table 6 and Table 4, respectively. To evaluate the accuracy of the constructed surrogate model, the coefficient of determination (R²) method was employed. A higher value, approaching 1, indicates stronger agreement between the predicted and actual values, and thus higher model accuracy. The calculation formula for R² is presented in Eq. (11).

R^{2} = 1 - \frac{\sum_{i = 1}^{n} {(y_{i} - {\hat{y}}_{i})}^{2}}{\sum_{i = 1}^{n} {(y_{i} - \bar{y})}^{2}}

(11)

RMSE = \sqrt{\frac{1}{n} \sum_{n}^{i = 1} {(F_{i} - {\hat{F}}_{i})}^{2}}

(12)

MAPE = \frac{100 %}{n} \sum_{n}^{i = 1} (F_{i} - {\hat{F}}_{i} / F_{i})

(13)

Table 6.

Simulated test values of composite laminates.

Test number	Stacking structure	Tensile strength (MPa)	Tensile modulus (GPa)	Elongation at peak (%)
1	CCCCCKCCC	542.28	180.66	5.68
2	CCCKCKCCC	518.61	172.56	5.55
3	CCCCCKKCK	541.10	175.93	6.08
4	CCKKCCCKC	460.97	165.81	5.95
5	CCKKKKCCC	476.11	159.23	5.28
6	KCKCCCKCK	455.35	157.77	5.12
7	CKKKCCKKC	424.71	155.91	5.76
8	CKKKKKKCC	465.99	154.96	6.41
9	KKCKKCKKC	431.38	143.28	5.86
10	KCKKKCKKK	455.53	144.56	6.47

where n is the number of sample points, $y_{i}$ is the true value of the $i$ -th sample, ${\hat{y}}_{i}$ is the Kriging predicted value of the $i$ -th sample, and $\bar{y}$ is the average of the true values of the samples.

Based on the results computed from Eq. (11), the R² values for tensile strength, tensile modulus, and elongation at break are 0.931, 0.919, and 0.917, respectively. In the context of engineering applications, an R² value exceeding 0.9 is generally considered sufficient for purposes. In addition, two additional error metrics, namely RMSE and MAPE, were adopted for a comprehensive evaluation of model prediction accuracy. The corresponding results are summarized in Table 7, which indicates a stable performance across all indicators. Therefore, the accuracy of the proposed surrogate model satisfies the criteria for subsequent optimization tasks (Wang et al., 2021). The distribution of the model fitting accuracy is illustrated in Fig. 7, which visually confirms the strong correlation between predicted and simulated tensile performance indicators.

Figure 7.

Comparison between predicted values and actual values (including experimental values and FEM values). (a) The tensile strength. (b) The tensile modulus. (c) The elongation at peak.

Table 7.

Prediction accuracy of three objective indicators.

Evaluation index	Tensile strength	Tensile modulus	Elongation at peak
R²	0.931	0.919	0.917
RMSE	11.15	3.54	0.11
MAPE (%)	2.01	1.89	1.44

Design and performance analysis of multi-objective optimization algorithms

Over the past decade, the continuous advancement of computational techniques has brought increased attention to heuristic and metaheuristic methods for solving complex optimization problems. Most of these methods are inspired by biological systems, physical laws, or mathematical principles (Mahmoudi et al., 2025). For example, PSO mimics the social behavior of flocking birds or schooling fish; the GA is inspired by biological evolution; and the SPEA is rooted in evolutionary strategies. These algorithms have been widely applied to tackle high-dimensional and nonlinear optimization problems. However, they have been primarily used for single-objective optimization, which limits their applicability in engineering problems involving multiple conflicting objectives. To overcome this limitation, multi-objective optimization algorithms have been developed to generate a set of trade-off (Pareto-optimal) solutions. In this study, a hybrid SPEA2–MOPSO algorithm is proposed, combining the advantages of SPEA2 and MOPSO. This hybrid algorithm is specifically designed to optimize the tensile properties of hybrid fiber-reinforced laminated composites.

Foundational algorithms and improvements

The PSO algorithm was originally developed by Kennedy and Eberhart (2001). To solve multi-objective problems, Coello et al. (2002) developed an improved version of the PSO algorithm, MOPSO, in 2004.

Algorithm 1 presents the pseudocode for MOPSO. The algorithm begins by initializing the particle swarm parameters and generating the initial positions and velocities of all particles in the population. It then evaluates each particle's objective functions and identifies the non-dominated solutions, which are used to construct an external archive that stores the current Pareto-optimal solution set. This archive is initially populated with the non-dominated solutions obtained from the initial swarm.

Algorithm 1.

MOPSO Pseudocode.

t \leftarrow 1

p \leftarrow

GenerateInitialSwarm()

A \leftarrow

Empty()

A \leftarrow

UpdateArchive(

p, A

)

5: while

t < t m a x

6: for each i in p do

x m g \leftarrow

SelectGlobalBest(

A

)

v (i) \leftarrow

UpdateVelocity(

x, x m p, x m g

)

x (i) \leftarrow

UpdatePosition(

x, v

)

10:

x (i) \leftarrow

Mutation(

x

)

11:

p (i) \leftarrow

CheckBoundary(

x

)

12: F

(i) \leftarrow

FitnessEvaluation(

x

)

13:

x m p (i) \leftarrow

UpdatePersonalBest(

F

)

14: end for

15:

A \leftarrow

UpdateArchive(

p, A

)

16:

t \leftarrow t + 1

17: end while

To clarify the operational framework of the algorithm, the following definitions are introduced. To maintain the diversity of the solution set and avoid premature convergence, the crowding distance is calculated and applied during archive truncation. The crowding distance is defined by Eq. (14) (Cao et al., 2021), where $f_{m} (i + 1)$ , $f_{m} (i - 1)$ , $f_{m}^{m a x}$ , and $f_{m}^{m i n}$ represent neighboring and boundary objective values.

d_{i} = \frac{f_{m} (i + 1) - f_{m} (i - 1)}{f_{m}^{m a x} - f_{m}^{m i n}}

(14)

Next, the position and velocity of each particle are updated according to Eqs. (15) and (16) (Shi and Eberhart 1998). A linearly decreasing inertia weight is incorporated to balance global and local search, defined in Eq. (17) (Nickabadi et al., 2011).

x_{i} (t + 1) = x_{i} (t) + ν_{i} (t)

(15)

ν_{i} (t + 1) = w ν_{i} (t) + c_{1} r_{1} (P_{i} (t) - x_{i} (t)) + c_{2} r_{2} (P_{b} (t) - x_{i} (t))

(16)

w = w_{\max} - (w_{\max} - w_{\min}) \frac{t}{T}

(17)

Here, $x_{i} (t)$ is the position vector of particle P at iteration t, and the velocity $ν_{i} (t)$ is added to the current position to update it. w is the inertia weight, the initial inertia damping coefficient is $W_{d a m p}$ = 0.99, $c_{1}$ and $c_{2}$ are the individual learning factor and the global learning factor, respectively. $r_{1}$ and $r_{2}$ are random numbers between 0 and 1. $P_{i}$ 、 $P_{b}$ represent the personal best and global best, respectively.

The original SPEA algorithm demonstrated notable performance in combinatorial optimization problems such as the 0/1 knapsack problem (Zitzler, 1999). SPEA2, proposed by Zitzler, Laumanns, and Thiele (Trentini et al., 2024), is an enhanced version incorporating density estimation and archive truncation to improve convergence and maintain diversity (Wu et al., 2024). Algorithm 2 provides the pseudocode for SPEA2.

SPEA2 introduces three key mechanisms that distinguish it from the original SPEA framework:

(1) Refined fitness assignment mechanism

The fitness of each individual in the population is computed based on two components. Strength value (S(i)): the number of other individuals dominated by individual i. Raw fitness(R(i)): It represents the number of other individuals that dominate individual i.

(2) Density estimation strategy

To enhance the distribution uniformity of solutions along the Pareto front, SPEA2 further incorporates local density estimation into the fitness assignment. The density estimator used in SPEA2 is based on the $k - t h$ nearest neighbor method, where the density of a point is determined by its distance to the $k - t h$ nearest solution. To determine the density D(i) of the $i - t h$ individual, it is necessary to calculate its distance to all other individuals j in both the archive and the population, forming a list $σ_{i}$ for each individual i that contains the distances to all other solutions j. During the distance evaluation process, each dimension in the objective space must be normalized according to the range of the objective space. After organizing $σ_{i}$ in ascending order, the $k - t h$ element of this list is denoted as $σ_{i}^{k}$ . Zitzler, Laumanns, and Thiele (Trentini et al., 2024) recommend using $k = \sqrt{N + \bar{N}}$ , where N is the population size and $\bar{N}$ is the archive size. Based on this, the density is defined as $D (i) = 1 / (σ_{i}^{k} + 2)$ . The final fitness is calculated as $F (i) = R (i) + D (i)$ .

(3) External archive management

SPEA2 employs a dual-population structure in which non-dominated solutions are stored in an external archive for use in the selection of the next generation. To preserve both the quality and diversity of the archived solutions, a density-based truncation strategy is applied when the number of non-dominated solutions exceeds the archive capacity. This strategy iteratively removes individuals located in densely populated regions, thereby preferentially retaining solutions situated in sparsely populated areas of the objective space.

Algorithm 2.

SPEA2 Pseudocode.

t \leftarrow 1

P \leftarrow

GenerateInitialPopulation()

F F \leftarrow

FineGrainedFitnessAassignment (

P

)

\bar{P} F \leftarrow

Survival(

P, F F

)

5: while

t < t m a x

F F \leftarrow

FineGrainedFitnessAassignment (

\bar{P}

)

M \leftarrow

MateSelection(

\bar{P}, F F

)

P \leftarrow

Crossover(

M

)

P \leftarrow

Mutation(

P

)

10:

P \leftarrow

BoundaryCheck(

P

)

11:

F \leftarrow

FitnessEvaluation (

P

)

12:

F F \leftarrow

FineGrainedFitnessAassignment (Union(

\bar{P,} P

))

13:

\bar{P} \leftarrow

Survival(Union(

\bar{P}, P), F F

)

14:

t \leftarrow t + 1

15: end while

During the environmental selection stage, SPEA2 first copies all non-dominated solutions (i.e., those with a strength fitness value less than 1) into the new archive for the next iteration. If the number of non-dominated solutions does not exceed the archive capacity, the environmental selection process terminates at this step. Otherwise, two situations may occur: the number of archived individuals may be smaller than the capacity, $| \bar{P_{t + 1}} | < \bar{N}$ , or it may exceed the capacity, $| \bar{P_{t + 1}} | > \bar{N}$ .

In the first case, the archive is supplemented by adding individuals with the smallest strength fitness values from the current population and the previous archive. In the second case, a truncation procedure is required to reduce the archive size to $\bar{N}$ .This truncation is performed iteratively: at each iteration, all pairwise distances are recalculated, and the individual with the smallest minimum distance is removed. In the case of a tie for the smallest distance in the $σ_{i}$ , the individual with the second smallest distance is removed, and so on.

This truncation mechanism ensures a well-distributed set of solutions across the objective space while guaranteeing that boundary solutions are never eliminated. The process of generating a new population, evaluating performance, and conducting environmental selection is repeated for a predefined number of generations. The final optimization results are represented by the solutions contained in the external archive.

Algorithm parameter settings

Parameter configuration constitutes a critical step in the optimization process, as parameter settings directly determine algorithmic performance and computational efficiency. Properly calibrated parameters can enhance solution quality, improve convergence stability, and reduce computational costs; conversely, inappropriate parameter selection may lead to inaccurate results or slow convergence. The importance of parameter selection has been validated across numerous engineering optimization problems. For example, Martínez-Muñoz et al., (2022) showed that parameter tuning in structural optimization applications can save up to 18% of computational time while achieving superior objective function values.

Many optimization algorithms, such as GAs, SPEA, and PSO, exhibit performance that is highly dependent on a set of internal parameters requiring calibration to achieve optimal results. For instance, within GA frameworks, parameters such as the crossover distribution index, mutation distribution index, and population size must be appropriately tuned. In the SPEA framework, convergence performance and solution-set distribution are primarily influenced by archive size, nearest-neighbor parameters used in fitness assignment, and the execution probabilities of genetic operators. Similarly, in PSO, the learning coefficients (C₁ and C₂) and the inertia weight (i) are key parameters governing its convergence behavior.

To systematically and efficiently determine the key parameter values for the hybrid algorithm employed in this study, the Taguchi method—a statistically grounded technique widely used to improve product and process quality in industrial engineering—was adopted. This method has also been successfully applied to the calibration of optimization algorithm parameters (Raz et al., 2021). Several studies have validated its effectiveness. For example, Mousavi et al. (2016) employed the Taguchi method to tune parameters of MOPSO and NSGA-II for inventory control optimization. Zeidabadi et al. (2023) used the method to calibrate SPEA2 and NSGA-II parameters in optimizing product delivery lead times. In addition, Trentini et al. (2024) conducted detailed parameter tuning for MOPSO, NSGA-II, and SPEA2 to solve engineering optimization problems, providing valuable insights for the present study. Based on these analyses and comparative evaluations, the final selection of key parameters adopted in this work is summarized in Table 8.

Table 8.

Parameters used for performance evaluation.

MOPSO
Parameter	Value	Parameter	Value	Parameter	Value
w	0.4	$W_{d a m p}$	0.99	$a l p h a$	0.1
C₁	0.75	$S_{w a r m}$	100	beta	2
C₂	1	$n_{G r i d}$	5	gamma	2
$η_{m}$	10
SPEA2
Parameter	Value	Parameter	Value	Parameter	Value
$η_{c}$	20	$η_{m}$	20	Population	100
$p_{c}$	0.8	$p_{m}$	1/9	$m_{u}$	0.1

Multi-objective optimization based on improved hybrid SPEA2-MOPSO algorithm

Given the complexity of the design variables involved in hybrid composite laminated plates, achieving both high optimization accuracy and computational efficiency is critical. To this end, the present study adopts the SPEA2 algorithm as the foundational framework and strategically integrates components of the MOPSO algorithm, dividing the overall optimization process into three distinct stages, as illustrated in Fig. 8. In the first stage (iterations < N/3), only the genetic operators of SPEA2—selection, crossover, and mutation—are employed. Mutation plays a key role in introducing diversity by expanding the search range, allowing for a broad exploration of the solution space and mitigating the risk of premature convergence to local optima. This phase focuses on rapidly constructing a high-quality, diverse set of initial non-dominated solutions, which serves as a strong foundation for subsequent search stages. In the second stage (N/3 < iterations<2N/3), the algorithm switches randomly between genetic operations and the velocity update mechanism of MOPSO. This hybrid strategy fuses the exploratory “jumping” behavior of genetic operators with the convergent, guided search characteristic of PSO. Non-dominated solutions stored in the archive act as leaders to guide the search, while random switching allows the algorithm to adapt dynamically to the current search state, thereby improving its robustness and adaptability. In the final stage (iterations>2N/3), the optimization relies solely on the MOPSO velocity update mechanism. At this point, the focus shifts to fine-tuning the Pareto front. The PSO mechanism drives particles toward both individual and global optima, while the inertia weight w decays linearly, reducing step size and thereby enhancing the local search precision. Non-dominated solutions in the archive continue to act as leaders, ensuring a uniform distribution across the Pareto front. This stage prioritizes solution refinement and convergence accuracy.

Figure 8.

The evolutionary methods at different stages of the hybrid SPEA2-MOPSO algorithm.

The initialization process of the hybrid algorithm is consistent with standard evolutionary approaches. It begins with the generation of the initial population and evaluation of corresponding response values using the previously developed surrogate model. The fitness value $F_{t}$ of each individual is then calculated using Eq. (18), which incorporates both objective function values and penalties for constraint violations:

F_{t} = \sum_{i = 1}^{2} f_{i} (x) + \sum_{j = 1}^{2} max (0, p_{j} (x))

(18)

where $f_{i} (x)$ denotes the objective functions (tensile strength, tensile modulus, elongation at break), $p_{j} (x)$ is the penalty function representing the degree of constraint violation.

After the initialization of the parent population, individual selection is performed using random tournament selection. Upon generating a new population, the MOPSO algorithm is introduced to guide search directions. Subsequently, the algorithm enters its three-stage hybrid optimization process, during which selection, crossover, and mutation operations are applied iteratively to generate the offspring population. The crossover operation is as shown in Eq. (19).

{\begin{matrix} y_{1} = 0.5 * [(1 + β) x_{1} + (1 - β) x_{2}] \\ y_{2} = 0.5 * [(1 - β) x_{1} + (1 + β) x_{2}] \end{matrix}

(19)

where

x_{1}

and

x_{2}

are the parent individuals,

y_{1}

and

y_{2}

are the offspring individuals, and

β

is the spread factor, which is determined by Eq. (20). u is a random number within the interval [0,1], and

η_{c}

is the crossover distribution index (

η_{c}

= 20), which controls the distribution range of offspring solutions around the parent solutions.

β = {\begin{matrix} {(2 u)}^{\frac{1}{η_{c} + 1}}, u \leq 0.5 \\ {[\frac{1}{2 (1 - u)}]}^{\frac{1}{η_{c} + 1}}, u > 0.5 \end{matrix}

(20)

To enhance individual diversity, polynomial mutation is applied to the individuals, with the specific formula given in Eq. (21). Here, $x_{m a x}$ and $x_{m i n}$ represent the boundary values of the variables, $δ$ denotes the mutation coefficient, and $η$ is set to 20.

\begin{matrix} x^{'} = x + δ (x_{max} - x_{min}) \\ δ = {\begin{matrix} {[2 u + (1 - 2 u) {(1 - δ_{1})}^{η + 1}]}^{\frac{1}{η + 1}} - 1, u ⩽ 0.5 \\ 1 - {[2 - 2 u + (2 u - 1) {(1 - δ_{2})}^{η + 1}]}^{\frac{1}{η + 1}}, u > 0.5 \end{matrix} \\ δ_{1} = (x - x_{\min}) / (x_{\max} - x_{\min}) \\ δ_{2} = (x - x_{\max}) / (x_{\max} - x_{\min}) \end{matrix}

(21)

The operations of merging the parent and offspring populations and generating a new population are repeated. Selection, crossover, and mutation operations are applied again to the new parent population to produce offspring populations until the maximum number of iterations is reached. The overall architecture and workflow of the proposed hybrid SPEA2–MOPSO algorithm is depicted in Fig. 9.

Figure 9.

Improved hybrid SPEA2-MOPSO algorithm flowchart.

The calculation of individual crowding distance follows Eq. (14). When the number of iterations is less than N/2, mutation is performed in accordance with Eq. (21). When the number of iterations exceeds N/2, mutation adheres to Eq. (22).

{\begin{matrix} Δ x_{N + 1}^{i} = Δ x_{N}^{i} + c_{1} r_{1} (x_p_{N}^{i} - x_{i}) + c_{2} r_{2} (x_g_{N} - x_{i}) \\ x_{i}^{'} = x_{i} + Δ x_{N + 1}^{i} \end{matrix}

(22)

where $x_g_{N}$ denotes the global optimal individual in the parent population at the $N - t h$ iteration, while $x_p_{N}^{i}$ represents the historical optimal individual at the $i - t h$ position in the $N - t h$ iteration. w stands for the inertia weight, $c_{1}$ denotes the individual learning weight, $c_{2}$ represents the population learning weight, and $r_{1}$ and $r_{2}$ are random numbers.

Comparative performance analysis of optimization algorithms

Although MOPSO can quickly converge by leveraging historical best solutions, it is prone to premature convergence and often struggles to maintain solution diversity, particularly in complex search spaces. SPEA2 shows strong global search ability and robustness in handling complex MOD problems with numerous discrete, non-dominated solutions. However, due to its large population and random evolutionary path, its convergence speed is relatively slow compared to PSO-based approaches.

Fig. 10 presents the optimization results obtained using the three algorithms. It can be observed that the optimized composite laminated plate exhibits notable improvements across all three tensile performance indicators—tensile strength, tensile modulus, and elongation at break. Among the compared algorithms, the SPEA2 algorithm produces solutions with relatively higher elongation at break, but lower tensile strength and modulus. This outcome suggests that SPEA2, while maintaining strong diversity, exhibits a slower convergence rate, particularly in optimizing tensile strength and modulus, making it difficult to reach ideal target values within a limited number of iterations.

Figure 10.

Pareto solutions under three algorithms.

By comparison, the hybrid SPEA2–MOPSO algorithm demonstrates superior overall performance. It successfully produces well-distributed and converged solutions, achieving a better balance among all three objectives. Consequently, the hybrid algorithm was selected for final implementation to optimize the tensile properties of the composite laminated plate. Fig. 11 illustrates the distribution of the Pareto-optimal solution set obtained by the hybrid algorithm. As shown in Fig. 11(a), most optimized solutions lie on the response surface, confirming the algorithm's ability to identify feasible, high-quality solutions. This not only highlights the hybrid algorithm effectiveness in multi-objective optimization, but also demonstrates the accuracy of the Kriging surrogate model employed in this study.

Figure 11.

The distribution of Pareto solution sets for the hybrid SPEA2-MOPSO algorithm. (a) 3D solution. (b) The tensile strength and tensile modulus. (c) The tensile modulus and elongation at peak. (d) The tensile strength and elongation at peak.

Further analysis of Fig. 11(b) reveals a strong positive correlation between tensile strength and tensile modulus, suggesting a near-linear relationship—an expected result due to the intrinsic mechanical interdependence of these properties. Fig. 11(c) and 11(d) show that the elongation at break values in the Pareto set remain relatively high across many solutions, reflecting the trade-off dynamics introduced during the optimization process. When higher elongation at break is prioritized, it inevitably results in a decrease in tensile strength and modulus, demonstrating the coupling effect among the three performance indicators. This interplay leads to an approximately linear inverse relationship between elongation at break and the other two tensile properties, further underscoring the complexity of balancing mechanical behaviors in laminated composite design.

Solution selection based on TOPSIS

Due to the inherent coupling effects among the three optimization objectives, the non-dominated solutions obtained represent a coordinated trade-off among competing performance indicators. For example, a relatively lower elongation at break may be acceptable if it accompanies higher tensile strength and modulus. Thus, the Pareto front offers a diverse set of feasible design alternatives for multi-objective problems. However, selecting a final solution from this set requires consideration of engineering requirements and application-specific constraints. Accordingly, this study discusses trade-off solutions in the final section.

To facilitate the decision-making process, the TOPSIS (Wang et al., 2016) is employed to transform the multi-objective optimization problem into a single-objective format. PCA (Xu and Wang, 2022) is subsequently incorporated to evaluate the Pareto-optimal solutions under varying weight distributions. The specific selection procedure is outlined as follows (Heidari et al., 2022):

Step 1. Construct the decision matrix A = $a_{i j}$ from the Pareto-optimal solution set.

Step 2. Normalize the matrix to obtain the normalized decision matrix B = $b_{i j}$ using Eq. (23).

b_{i j} = a_{i j} / \sqrt{\sum_{i = 1}^{m} a_{i j}^{2}}

(23)

Step 3. Calculate the weighted decision matrix C = $c_{i j}$ as shown in Eq. (24).

c_{i j} = w_{j} \cdot b_{i j}

(24)

where

w_{j}

is the weight of the j objective.

Step 4. Determine the positive ideal solution (PIS) and negative ideal solution (NIS) based on Eqs. (25) and (26).

P I S = max_{i} c_{i j}

(25)

N I S = min_{i} c_{i j}

(26)

Step 5. Compute the Euclidean distance of each non-dominated solution from the PIS and NIS using Eqs. (27) and (28).

d_{i}^{P} = \sqrt{\sum_{j = 1}^{n} {(c_{i j} - P I S)}^{2}}

(27)

d_{i}^{N} = \sqrt{\sum_{j = 1}^{n} {(c_{i j} - N I S)}^{2}}

(28)

Step 6. Calculate the closeness coefficient $R_{i}$ for each solution using Eq. (29), and rank the solutions accordingly. The closer $R_{i}$ is to 1, the higher the ranking, and the solution is selected as the optimal candidate.

R_{i} = d_{i}^{N} / (d_{i}^{N} + d_{i}^{P})

(29)

After transforming the multi-objective optimization problem into a single-objective format, the optimal solution is selected by analyzing the influence of each weight coefficient on the variation of the three tensile performance indicators under different weight distributions. Fig. 12 illustrates the sensitivity of these indicators to different weight allocations, where $w_{1}$ and $w_{2}$ denote the weights assigned to tensile strength and elongation at break, respectively. As shown in Fig. 12(a), when $w_{2}$ is small, tensile strength is highly responsive to changes in $w_{1}$ . Specifically, as $w_{1}$ increases, the tensile strength steadily improves, indicating that $w_{1}$ is the dominant factor influencing this property. Conversely, when $w_{1}$ is small, increasing $w_{2}$ initially leads to a slight improvement in tensile strength, which then plateaus, suggesting that the effect of $w_{2}$ on tensile strength is relatively limited. Fig. 12(b) illustrates the effect of weight coefficients on tensile modulus. When either $w_{1}$ or $w_{2}$ is small, the tensile modulus increases with the increase of the other coefficient. However, when either weight becomes large, the tensile modulus shows little variation. This indicates that the tensile modulus is primarily governed by its own corresponding weight coefficient rather than by $w_{1}$ or $w_{2}$ . Similarly, Fig. 12(c) presents the weight sensitivity analysis for elongation at break. When $w_{1}$ is small, increasing $w_{2}$ leads to a gradual increase in elongation at break, confirming the strong influence of $w_{2}$ on this parameter. In contrast, when $w_{2}$ is small, changes in $w_{1}$ have minimal effect, further reinforcing that elongation at break is more sensitive to its own weight coefficient than to that of tensile strength.

Figure 12.

The effect of weight coefficient. (a) Tensile strength. (b) Tensile modulus. (c) Elongation at peak.

Due to variations in engineering design requirements, the final decision must be made based on specific application scenarios. Table 9 summarizes the selection results for a subset of cases. For the first three solutions listed, tensile strength, tensile modulus, and elongation at break were designated as the primary optimization objectives, respectively. In each case, the main objective was assigned a weight of 0.6, while the remaining two objectives were weighted at 0.2. In contrast, Solution 4 considered all three objectives equally important, assigning equal weights to each. The distribution of selected solutions from the Pareto front is illustrated in Fig. 13.

Figure 13.

The Pareto optimal selection.

Table 9.

Pareto optimal feasible solutions for different decisions.

Solutionnumber	Stackingstructure	Tensile strength(MPa)	Tensile modulus(GPa)	Elongation at peak(%)
1	CCCCCCCKK	577.35	194.23	6.01
2	CKCCCCCCK	586.87	197.41	6.02
3	CKCCKKKCK	547.43	185.10	6.72
4	KCKCCCCKC	572.25	192.48	6.27

The hybrid SPEA2-MOPSO algorithm significantly improved the comprehensive tensile performance of the optimized composite laminates. Compared with pure CFRP laminates, the optimized solutions exhibit reduced tensile strength but enhanced elongation at break. Conversely, compared with pure KFRP laminates, the optimized solutions show substantial improvements in tensile strength, accompanied by a moderate reduction in elongation at break. This trend is attributed to the distinct characteristics of the base materials: pure CFRP exhibits high tensile strength but limited ductility, while pure KFRP offers excellent ductility but relatively low strength. The optimized laminates, however, achieve a more balanced performance across all three indicators, namely tensile strength, modulus, and elongation, demonstrating the effectiveness of the proposed hybrid optimization algorithm.

Compared with the CF5/KF4 laminate before optimization, in Solution 1, the elongation at break decreased by 6.2%, while tensile strength and modulus increased by 7.6% and 9.1%, respectively. Solutions 2 and 4, which emphasized strength and modulus, achieved even higher improvements in those metrics, with relatively minor reductions in elongation at break (6.1% and 2.2%, respectively). In contrast, solution 3 achieved simultaneous improvements in all three performance metrics, namely tensile strength, modulus, and elongation at break, increasing by 2.0%, 3.9%, and 4.8%, respectively. In summary, although the optimized solutions differ in the degree of improvement across individual performance metrics, all show enhanced comprehensive tensile properties, validating the correctness and effectiveness of the proposed optimization framework. It is important to note that the four feasible solutions presented here serve as illustrative examples to demonstrate the advantages of the hybrid SPEA2-MOPSO algorithm optimization strategy compared with the baseline laminate. In practical engineering applications, the weighting of each objective should be adjusted based on actual design priorities to meet specific functional requirements.

Conclusion

This study presents a practical, efficient, and scalable methodology for the structural optimization of hybrid composite laminates. The validated high-fidelity simulation model, combined with the Kriging surrogate model and hybrid SPEA2-MOPSO algorithm, offers a robust framework for exploring design strategies under multi-criteria constraints. Additionally, the integration of the TOPSIS method enhances decision-making by enabling tailored selection based on application-specific priorities. The main conclusions are summarized as follows:

Tensile experiments and progressive damage finite element simulations were conducted on carbon/Kevlar hybrid laminates to validate the accuracy of the developed model. Simulations and tests across various stacking sequences and material combinations yielded maximum relative errors within 7.18% for tensile strength, tensile modulus, and elongation at break. The stress–strain curves from simulations closely matched experimental data, with a maximum relative error of 6.48% in tensile modulus, confirming the high fidelity of the finite element model.

The optimal super LHS method, in combination with the Kriging surrogate model, enabled accurate prediction of laminate tensile performance. A database of design samples was generated via numerical simulation, capturing tensile property metrics and reliability parameters. The determination coefficient (R²) values for the test set exceeded 0.9 across all three performance indicators, demonstrating both the uniformity of the sampling approach and the predictive accuracy of the surrogate model.

The hybrid SPEA2-MOPSO algorithm effectively addressed the multi-objective optimization problem for composite laminate tensile properties. Compared with standalone SPEA2 and MOPSO algorithms, the hybrid method demonstrates superior performance and optimization efficiency, offering key advantages in improving predictive accuracy and reducing computational cost. Furthermore, the algorithm generated diverse Pareto-optimal solution sets, enhancing both robustness and applicability in engineering design scenarios by offering broader selection flexibility.

To address varying engineering requirements, the TOPSIS, combined with PCA, was employed to assign weights and select the most suitable compromise solutions from the Pareto front. The results revealed a well-balanced improvement across tensile strength, tensile modulus, and elongation at break, confirming the effectiveness of the proposed optimization strategy.

Overall, the proposed framework not only demonstrates significant potential in improving the tensile performance of composite laminates but also provides a valuable reference for tackling more complex structural optimization problems. The development and application of such multi-objective optimization techniques are pivotal for advancing composite material design and broader structural engineering applications.

Footnotes

Acknowledgments

The authors are grateful for support from the Natural Science Foundation Project of the Fujian Provincial Department of Science and Technology (Grant No. 2024J01130754), the Natural Science Foundation Project of the Xiamen Municipal Science and Technology Bureau (Grant No. 3502Z202471076). This work was supported by open fund of Fujian Provincial Key Laboratory of Advanced Design and Manufacture for Bus Coach (Xiamen University of Technology). The authors also sincerely extend their gratitude to the editors and reviewers for providing revision suggestions.

ORCID iD

Jiayou Guo

CRediT authorship contribution statement

Hangyan Wang: Experiment, Investigation, Writing-original draft, Writing-review & Editing. Jiayou Guo: Conceptualization, Data curation, Writing-original draft & Investigation. Hao Fang and Junjie Lin: Investigation, Project administration & Supervision.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Natural Science Foundation Project of the Xiamen Municipal Science and Technology Bureau, Natural Science Foundation Project of the Fujian Provincial Department of Science and Technology, (grant number Grant No. 3502Z202471076, Grant No. 2024J01130754).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data availability statement

Data will be made available upon reasonable request.

References

Arhore

Yasaee

(2020) Lay-up optimisation of fibre–metal laminates panels for maximum impact absorption. Journal of Composite Materials 54(29): 4591–4609.

ASTM D3039 . (n.d.) Standard tests method for tensile properties of polymer matrix composite materials.

Azhdari

Taheri-Behrooz

. (2023). An experimental and numerical investigation on bending after impact strength of glass laminate aluminium reinforced epoxy. Composites Part A-Applied Science And Manufacturing 2023: 171.

Berton

Coussa

Berthe

, et al. (2026) Investigation of strain rate effects on the irreversible and damageable behavior of carbon fiber reinforced polymers under dynamic loading. International Journal of Damage Mechanics 35(1): 3–21.

Beylergil

(2020) Multi-objective optimal design of hybrid composite laminates under eccentric loading. Alexandria Engineering Journal 59(6): 4969–4983.

Bresciani

Manes

Ruggiero

, et al. (2016) Experimental tests and numerical modelling of ballistic impacts against kevlar 29 plain-woven fabrics with an epoxy matrix: Macro-homogeneous and meso-heterogeneous approaches. Composites Part B: Engineering 88: 114–130.

Cao

Khan

Abdi

, et al. (2021) Combination of RSM and NSGA-II algorithm for optimization and prediction of thermal conductivity and viscosity of bioglycol/water mixture containing SiO2 nanoparticles. Arabian Journal of Chemistry 14(7): 103204.

Chairi

El Bahaoui

Hanafi

, et al. (2026) Computational assessment of sustainable balsa and basalt composite sandwich for structural marine applications. International Journal of Damage Mechanics 35(2): 177–200.

Chen

Bhandari

, et al. (2023) Genetic optimized al–mg alloy constitutive modeling and activation energy analysis. International Journal of Mechanical Sciences 244: 108077.

10.

Cheng

Yuan

, et al. (2020) Transverse and longitudinal flexural properties of unidirectional carbon fiber composites interleaved with hierarchical aramid pulp micro/nano-fibers. Composites Part B: Engineering 188: 107897.

11.

Coello

Ing

Lechuga

. (2002) MOPSO: A Proposal for Multiple Objective Particle Swarm. In Proceedings of the 2002 Congress on Evolutionary Computation, IEEE 2002; CEC'02 (Cat. No. 02TH8600) (Vol. 2, pp. 1051-1056).

12.

Fakoor

Ghoreishi

SMN

Aminjafari

(2019) Multi-objective optimization of buckling load for a laminated composite plate by coupling genetic algorithm and FEM. Journal of Aerospace Science and Technology 12(1): 27–37.

13.

Garg

Mukhopadhyay

Belarbi

, et al. (2023) On accurately capturing the through-thickness variation of transverse shear and normal stresses for composite beams using FSDT coupled with GPR. Composite Structures 305: 116551.

14.

Garg

Mukhopadhyay

Belarbi

, et al. (2023) Random forest-based surrogates for transforming the behavioral predictions of laminated composite plates and shells from FSDT to elasticity solutions. Composite Structures 309: 116756.

15.

Ghane

Fagerström

Mirkhalaf

(2024) Recurrent neural networks and transfer learning for predicting elasto-plasticity in woven composites. European Journal of Mechanics-A/Solids 107: 105378.

16.

Ghiasi

Pasini

Lessard

(2009) Optimum stacking sequence design of composite materials part I: Constant stiffness design. Composite Structures 90(1): 1–11.

17.

Gholami

Fathi

Baghestani

(2021) Multi-objective optimal structural design of composite superstructure using a novel MONMPSO algorithm. International Journal of Mechanical Sciences 193: 106149.

18.

Gower

Cronin

Plumtree

(2008) Ballistic impact response of laminated composite panels. International Journal of Impact Engineering 35(9): 1000–1008.

19.

Heidari

Boleydei

Rohani

, et al. (2022) Integrating life cycle assessment and life cycle costing using TOPSIS to select sustainable biomass-based-carbonaceous adsorbents for CO2 capture. Journal of Cleaner Production 357: 131968.

20.

Hoshikawa

Ryuzono

Onodera

, et al. (2025) Finite element modeling of viscoelastic creep behavior and transverse cracking in fiber-reinforced composite materials. International Journal of Damage Mechanics 34(9): 1365–1387.

21.

Hou

Zhang

, et al. (2023) Design and optimization of the bumper beam with corrugated core structure of fiber metal laminate subjected to low-velocity impact. Thin-Walled Structures 187: 110746.

22.

Houari

Siguerdjidjene

Chellil

, et al. (2026) Experimental and numerical study of the vibration behavior of CFRP plates repaired with a bonded composite patch and its adhesive. Journal of Vibration Engineering & Technologies 14: 64.

23.

Huang

Tao

Sun

, et al. (2023) Assessment of numerical modeling approaches for thin composite laminates under low-velocity impact. Thin-Walled Structures 191: 111053.

24.

Imran

Shi

Tong

, et al. (2021) Multi-objective design optimization of composite submerged cylindrical pressure hull for minimum buoyancy factor and maximum buckling load capacity. Defence Technology 17(4): 1190–1206.

25.

Jiang

Ren

(2022) CFRP-patching enhancement on open-hole CFRP panel with micro/nanofillers-modified adhesive interface: Experimental and numerical simulation. Composites Science and Technology 218: 109180.

26.

Jiang

Cheng

, et al. (2021) Micro-mechanics modeling of compressive strength and elastic modulus enhancements in unidirectional CFRP with aramid pulp micro/nano-fiber interlays. Composites Science and Technology 206: 108664.

27.

Jiang

Ren

Jin

, et al. (2020) Flexural performances of fiber face-sheets/corrugated core sandwich composite structures reinforced by horizontal stiffeners. International Journal of Mechanical Sciences 168: 105307.

28.

Kennedy

Eberhart

(2001) Swarm Intelligence. Burlington, MA: Morgan Kaufmann Publishers.

29.

Liu

, et al. (2019) Assessment of failure criteria and damage evolution methods for composite laminates under low-velocity impact. Composite Structures 207: 727–739.

30.

, et al. (2023) Buckling analysis of variable-stiffness composite truncated elliptical rotary shell under external hydrostatic compression. Composite Structures 307: 116613.

31.

Mahmoudi

Maleki

Rezaei Ochbelagh

(2025) Multi-objective optimization of hybrid energy systems using gravitational search algorithm. Scientific Reports 15(1): 2550.

32.

Martínez-Muñoz

García

Martí

, et al. (2022) Optimal design of steel–concrete composite bridge based on a transfer function discrete swarm intelligence algorithm. Structural and Multidisciplinary Optimization 65: 12.

33.

Mir

MSS

Rezaeian

(2016) A robust hybrid approach based on particle swarm optimization and genetic algorithm to minimize the total machine load on unrelated parallel machines. Applied Soft Computing 41: 488–504.

34.

Mostofi

Sayah-Badkhor

Rezasefat

, et al. (2020) Gas mixture detonation load on polyurea-coated aluminum plates. Thin-Walled Structures 155: 106851.

35.

Mousavi

Sadeghi

Niaki

STA

, et al. (2016) A bi-objective inventory optimization model under inflation and discount using tuned pareto-based algorithms: NSGA-II, NRGA, and MOPSO. Applied Soft Computing 43: 57–72.

36.

Nickabadi

Ebadzadeh

Safabakhsh

(2011) A novel particle swarm optimization algorithm with adaptive inertia weight. Applied Soft Computing 11(4): 3658–3670.

37.

Peng

Chen

Jia

, et al. (2025) Stacking sequence optimization of variable thickness composite laminated plate based on multi-peak stacking sequence table. Composite Structures 356: 118886.

38.

Raz

Akabarzadeh

TMR

Akabarzadeh

(2021) Experiment-based affect heuristic using fuzzy rules and Taguchi statistical method for tuning complex systems. Expert Systems with Applications 172: 114638.

39.

Sahu

Andhare

(2019) Multiobjective optimization for improving machinability of Ti-6Al-4 V using RSM and advanced algorithms. Journal of Computational Design and Engineering 6(1): 1–12.

40.

Shi

Eberhart

(1998) A modified particle swarm optimizer. In 1998 IEEE international conference on evolutionary computation proceedings. IEEE world congress on computational intelligence (Cat. No. 98TH8360); 1998, May ;(pp. 69-73).

41.

Slamani

Chatelain

(2015) Issues and challenges in robotic trimming of CFRP. In: 2015 12th International Conference on Informatics in Control, Automation and Robotics (ICINCO). IEEE, 400–405.

42.

Sokolinsky

(2011) Numerical simulation of the crushing process of a corrugated composite plate. Composites Part A-Applied Science And Manufacturing 42(9): 1119–1126.

43.

Trentini

EVW

Parsekian

Bittencourt

(2024) Multiobjective optimization of bridge and viaduct design: Comparative study of metaheuristics and parameter calibration. Engineering Structures 312: 118252.

44.

Tsai

Liu

Chou

(2004) Hybrid taguchi-genetic algorithm for global numerical optimization. IEEE Transactions on Evolutionary Computation 8(4): 365–377.

45.

Wang

Cai

(2018) Multi-objective crashworthiness optimization of vehicle body using particle swarm algorithm coupled with bacterial foraging algorithm. Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering 232(8): 1003–1018.

46.

Wang

Jiang

(2016) A hybrid method of modified NSGA-II and TOPSIS for lightweight design of parameterized passenger car sub-frame. Journal of Mechanical Science and Technology 30: 4909–4917.

47.

Wang

Guo

Liu

(2025) Tensile properties and failure of hybrid fiber reinforced polymer composite laminate with different fiber types, hybrid ratios, and stacking sequences. Polymer Composites 46(10): 9348–9364.

48.

Wang

Guo

Zhang

, et al. (2024) Dynamic responses and interactive failure mechanisms of carbon fiber composite face sheets/double-layer corrugated core sandwich structures under low-velocity impacts loading. Journal of Composite Materials 58(14): 1695–1713.

49.

Wang

Lin

Jiang

, et al. (2024) Inter-layer failure and toughening mechanisms of carbon/aramid hybrid fiber composites interleaved with micro/nano pulps under low-velocity impact load. Thin-Walled Structures 202: 112086.

50.

Wang

Liu

Ren

, et al. (2023) Hybrid effects and failure mechanisms of carbon/kevlar fibre composite laminates under the bending-after-impact loading. Engineering Structures 299: 117101.

51.

Wang

Zhang

Zhou

, et al. (2021) Implementation of a novel six sigma multi-objective robustness optimization method based on the improved response surface model for bumper system design. Thin-Walled Structures 167: 108257.

52.

Wang

Pan

Ying

, et al. (2022) Symmetric and asymmetric intercalation effect on the low-velocity impact behavior of carbon/Kevlar hybrid woven laminates. Composite Structures 297: 115919.

53.

Liu

Yan

, et al. (2024) Blast performance of polyurethane foam-filled auxetic honeycomb sandwich beams. Composite Structures 338: 118104.

54.

Wang

(2022) Fatigue/impact analysis and structure-connection-performance integration multi-objective optimization of a bolted carbon fiber reinforced polymer/aluminum assembled wheel. Composites Part B: Engineering 243: 110103.

55.

Chen

Cui

, et al. (2024) Low-velocity impact analysis and multi-objective optimization of hybrid carbon/basalt fibre reinforced composite laminate. Composite Structures 343: 118305.

56.

Zeidabadi

Nik

Hejazi

(2023) Multi-objective robust production planning considering workforce efficiency with a metaheuristic solution approach. International Journal of Industrial Engineering 30(1): 32–50.

57.

Zhang

Chen

Zhuo

, et al. (2025) Machine learning-based optimization of impact-resistant layups for FRP hybrid laminates in yachts. Ocean Engineering 324: 120600.

58.

Zhang

Wang

Song

, et al. (2025) Ballistic performance and thickness optimization of TC4/UHMWPE laminates composite structure under high-speed impact. Composites Communications 56: 102411.

59.

Zhang

Zhou

, et al. (2022) Space deployable bistable composite structures with C-cross section based on machine learning and multi-objective optimization. Composite Structures 297: 115983.

60.

Zhao

Zhu

Zhou

, et al. (2019) Crashworthiness analysis and design of composite tapered tubes under multiple load cases. Composite Structures 222: 110920.

61.

Zheng

Sun

, et al. (2017) Tensile response of carbon-aramid hybrid 3D braided composites. Materials & Design 116: 246–252.

62.

Zhou

Wen

Wang

(2020) Numerical investigation on the repeated low-velocity impact behavior of composite laminates. Composites Part B: Engineering 185: 107771.

63.

Zitzler

(1999) Evolutionary Algorithms for Multiobjective Optimization: Methods and Applications. Shaker Ithaca.