A machine learning model for predicting the ballistic impact resistance of unidirectional fiber-reinforced composite plate

It has been a vital issue to ensure both the accuracy and efficiency of computational models for analyzing the ballistic impact response of fiber-reinforced composite plates (FRCP). In this paper, a machine learning (ML) model is established in an effort to bridge the ballistic impact protective performance and the characteristics of microstructure for unidirectional FRCP (UD-FRCP), where the microstructure of the UD-FRCP is characterized by the two-point correlation function. The results showed that the ML model, after trained by 175 cases, could reasonably predict the ballistic impact energy absorption of the UD-FRCP with a maximum error of 13%, indicating that the model can ensure both computational accuracy and efficiency. Besides, the model’s critical parameter sensitivities are investigated, and three typical ML algorithms are analyzed, showing that the gradient boosting regression algorithm has the highest accuracy among these algorithms for the ballistic impact problem of UD-FRCP. The study proposes an effective solution for the traditional difficulty of the ballistic impact simulation of composites with both high efficiency and accuracy.

Fiber-reinforced composite plate (FRCP) has been widely used in engineering involving aircraft structure, civil construction, ocean engineering, etc., due to its high specific strength and modulus [1][2][3][4][5] . In addition, the numerous energy dissipation channels including tension and failure of fibers, friction between fibers and matrix, and compression failure of matrix render FRCP excellent impact energy dissipation capacity [6][7][8] , making it a preferred material in impact protective structures such as bulletproof vests and vehicles.
The impact resistance of FRCP is dependent on the dynamic behavior of fibers and matrix, and it is affected significantly by the microstructure. Generally, there is a consensus that unidirectional FRCP (UD-FRCP) can provide higher impact resistance than other structural topologies 9 for the bulletproof composite. In addition to ballistic impact tests, numerical simulation has been proven to be an effective method for predicting the impact resistance of UD-FRCP [10][11][12] . The simulation method can be classified into two scales, i.e. macro-scale 13 and micro-scale 14,15 , according to the modeling characteristics. For the macro-scale simulation of FRCP implemented by the finite element analysis (FEA), the transverse isotropic constitutive model is usually used for each layer of UD-FRCP in cooperation with some failure criteria such as Hashin failure 16 to analyze the ballistic impact responses. For this method, the results such as ballistic limit velocity and damage morphology always agree with the experimental results according to the study by Abdel-Nasser et al. 17 , Zhu et al. 18 , etc. However, it is hard to obtain the evolution of composite microstructures during impact like damage behavior induced by the ballistic impact. To capture the detailed deformation and failure behavior of composites to increase the accuracy of computation, the micro-scale simulations are generally performed, in which fiber bundles, matrix, and architectures are modeled. According to Refs. 19,20 that are simulated by micro-scale models, the failure criteria of the matrix and the interfacial condition play a dominant role in the phenomenological failure criteria of the composite. The failure behaviors and the damage evolution of the composites are also clearly delineated during penetration. However, it is worth noting that the computation at micro-scale generally takes more than 2000 h CPU-time for one case 21 , which is extremely time-consuming when compared to the macro-scale simulation. As can be seen from the literature mentioned above, the computational efficiency of the macro-scale models is high, whereas the accuracy is generally not satisfied. In contrast, the computational accuracy of the micro-scale models is usually www.nature.com/scientificreports/ satisfactory, but the computational cost is too expensive for a model with relatively large dimensions. Therefore, it is crucial to put forward computation models with both high efficiency and accuracy for predicting the ballistic impact resistance of UD-FRCP.
In recent years, with the rapid development of computing science and technology, machine learning (ML) method begin to be widely adopted to predict the mechanical behavior of various materials, providing a promising solution for obtaining a computation model with both high efficiency and accuracy to obtain the ballistic impact resistance of UD-FRCP. The availability of supervised ML and convolutional neural networks have been verified in depicting the relationship between materials, structures, and mechanical properties 22,23 such as stress-strain curve, modulus of elasticity, and yield strength. The development of ML greatly broadens the research method of mechanics. However, the previous prediction studies by the ML method mainly focus on the quasi-static mechanical properties of composites. Whether the ML method can reasonably predict the ballistic impact behavior of UD-FRCP is still unknown.
In this paper, the ballistic impact resistance of UD-FRCP is studied by a ML method for the first time. The critical parameter sensitivities, various ML algorithms, and the influence of the number of training cases are also investigated. The results show that the ballistic impact resistance of UD-FRCP can be well predicted by the proposed model. The paper is organized as follows. "Model and methodology" gives the details of the ML model. In "Results and discussion", the prediction results are provided. The sensitivities of parameters are investigated, the effects of various ML algorithms are also compared, and the influence of the number of cases on the prediction error is demonstrated, followed by the discussion and conclusions.

Model and methodology
The detailed process of the ML model is shown in Fig. 1. Firstly, various microstructures of UD-FRCP are generated and characterized. The ballistic impact database of the UD-FRCP with different microstructure is then built based on massive FEA simulations. After that, three typical ML algorithms are employed to predict the ballistic impact resistance of the UD-FRCP with a given microstructure characteristic, and the sensitivities of the related parameters are investigated. Finally, a ML model with both high efficiency and accuracy is proposed for predicting the ballistic impact resistance of UD-FRCP.
Numerical model. Construction of UD-FRCP. We randomly generated 185 UD-FRCP materials with different cross-section characteristics with 18 ~ 48% fiber content. The method of microstructure generation is shown in Fig. 2a. Several microstructures are generated for the same fiber content in order to explore the influence of the fiber distribution characteristics on the impact protective performance of the composite materials.  Micro-scale ballistic impact model. The ABAQUS explicit software, which has been widely used in the impact response of composites [25][26][27] , is employed to analyze the ballistic impact response of the UD-FRCP. As shown in Fig. 3a, a 1/4 model with the dimensions of 100 × 100 × 10 mm is built due to the symmetries of the problem to improve computational efficiency. In addition, since hundreds of impact processes will be simulated to create the database for ML, a 1/8 rigid spherical surface with a mass of 55.8 g and a radius of 15 mm is taken as a projectile to mimic the bullet to decrease the computation cost further. The mesh sizes of the fibers and the matrix are matched, which are both 0.16 mm in the cross-sectional direction. In the extensional direction of the length, the size gradually transforms from 0.5 mm to 6.65 mm in order to reduce the number of meshes and consider the main effects of numerically simulated impact position. The meshes of the fiber bundle and the matrix exceed 620,000 and 640,000 C3D8R elements in the 1/4 model, respectively. The initial impact velocity of the projectile is 300 m/s. The symmetrical boundary conditions are applied to the symmetrical cross-sections of the UD-FRCP and the projectile, and the fixed boundary condition is applied to the peripheries of the UD-FRCP far away from the projectile. The general contact condition is used to model the non-linear contact behavior between the projectile and the UD-FRCP during penetration. Honestly, the ballistic impact models for the UD-FRCP are significantly simplified in an effort to achieve the highest computational efficiency. However, the effects of microstructures of the cross-sections on the impact protective performance of the UD-FRCP can still be revealed. One could expect better results from fine-meshed practical micro-scale models. Nevertheless, it will be quite time-consuming while building the database for ML, which will be performed in the future. Here, we mainly focus on the validation of the ML method for analyzing the impact protective performance of UD-FRCP with various microstructures.
The Tsai-Wu failure criterion 28 is undertaken in ABAQUS explicit package 29 to describe the failure behavior of the UD-FRCP, The failure parameters are defined as where X and X', Y and Y', Z and Z' are tensile and compressive strengths of the UD-FRCP along the fiber axial direction (x 1 ), the transverse direction (x 2 ), and the thickness direction (x 3 ), respectively; S 12 , S 13 , and S 23 are the shear strengths of the UD-FRCP in the x 3 , x 2 , and x 1 planes, respectively. The related material parameters of the fiber and the matrix are given in Table 1 Ref. 30 .
(2) www.nature.com/scientificreports/ The typical deformation and failure behavior of the UD-FRCP after penetration is shown in Fig. 3b, where the full model is re-constructed during post-processing to show the deformation and failure behavior of the UD-FRCP clearly. The residual impact velocities of the projectiles after penetrating the UD-FRCP with various cross-section characteristics can be obtained. Therefore, the dissipated impact energy E , which is used to evaluate the impact protective performance of the UD-FRCP, can be calculated as where m, v 0, and v r are mass, initial velocity, and residual velocity of the projectile, respectively.
For each impact process simulation, it takes nearly 2 h for a computer with an Intel (R) Core (TM) i7-8700 CPU and a memory of 16 GB. As a result, 185 impact processes of the UD-FRCP with various microstructures are simulated and analyzed to build the training set and the testing set to balance the ML accuracy and the computational cost. More cases will be beneficial to improve ML accuracy, but more time-consuming.

Microstructure characterization. Two-point correlation function.
It is necessary to extract the specific characteristic data of the cross-sections' microstructures from the cross-section images of the UD-FRCP, as given in Fig. 2b. Here, the n-point spatial correlation function, which is a measurement of the probability of finding n points in the space occupied by a component within a two-phase material, is employed to analyze the microstructure characteristics of the UD-FRCP. In detail, the one-point correlation function reveals the probability of any point in a type of material; the two-point correlation function reflects the likelihood of two points in a certain distance at the same time; the three-point correlation function indicates the probability of a specific triangle in a particular material. For the present problem, the two-point correlation function is applicable to reveal the microstructure characteristics of the UD-FRCP, such as the fiber content and the relative positions of fiber bundles in the cross-section 31 . The two-point correlation function defines the probability of finding phase P and phase P' at the head and the tail, respectively, of a vector with a length of r that is randomly placed in a microstructure represented by a voxel image with a × b pixels, as shown in Fig. 4. Each pixel in the image is identified by a unique two-dimensional position vector s. Therefore, the two-point correlation function is defined as where m P s is microstructure function taking the value of either 1 if the material phase P is located in the spatial position s or 0 if any other different phase is present at the spatial position s. S denotes the total number of pixels in the image 22 .
In this paper, the inter-fiber spatial correlation is analyzed, where F denotes fiber phase.
The microstructure of the composite in the present study is analyzed by 100,000 (1000 × 100) pixels. Figure 5a shows five examples of microstructures with 100, 115, 130, 130, and 130 fiber bundles corresponding to the fiber contents of 31.41%, 36.11%, 40.84%, 40.84%, and 40.84%, respectively. The corresponding two-point correlation  www.nature.com/scientificreports/ matrices (1000 × 100 pixels) of the microstructures are calculated as depicted in Fig. 5b. For the vector r with a length of 0 pixel, the value of f FF r represents the content of fibers. The two-point correlation function matrix differs significantly for different fiber contents. Even for the same fiber content as given in Fig. 5, the twopoint correlation functions show obvious differences due to the variation of the microstructures. Therefore, the two-point correlation function matrix can describe the characteristics of the composite microstructure. Each microstructure consists of 100,000 values, and a substantial two-correlation function matrix with a dimension of 100,000 × 185 is finally constructed for the total 185 cases.
Dimensionality reduction. Since the two-correlation function matrix is massive, it is tough to train the ML model based on the data of the matrix. Therefore, it is essential to reduce the dimensionality but retain the characteristics of the data set. Here, the principal component analysis (PCA) 32 , which is helpful for extracting the main features and reducing the size of the data set, is taken to shrink the dimensionality of the data set of the two-correlation function matrix. The principal components of a data set are found through linearly transform from a series of possible related variables to some unrelated variables. To ensure that the first principal component describes the direction of maximum variance, the data set has been centralized. According to the zero empirical mean theory 33 , the first weight vector w 1 of the data set X can be defined as which represents the maximum variance of w ⊤ X when taking w = w 1 .
For finding the kth principal component, where k ≥ 2, the previously obtained principal components need to be removed from the data set X, and this process can be depicted as follow, The minimum k is determined according to the maximum mean-variance of the data set, where m is the total data of the data set, x (i) and x (i) approx are the ith actual value and estimated value. e can be chosen by specific demand according to the trade-off between dimension and precision. The smaller the dimension is, the less accuracy in model training is achieved. In the present study, e is taken as 0.01 to ensure adequate accuracy for the main characteristics of the microstructures, and therefore k is determined as 10.
Machine learning algorithms. Gradient boosting regression algorithm. The gradient boosting regression (GBR) algorithm combines "weak learners" into a "strong learner" in an iterative manner 34 . At each stage of gradient boosting m, 1 ≤ m ≤ M, assuming there is already an imperfect model F m , the next-step model F m+1 of the GBR algorithm is improved by adding a new estimator h 35 , Therefore, gradient lifting obtains h by fitting y-F m (x). Like other variants of lifting methods, F m+1 becomes more accurate by correcting the error of F m . The idea is originated from the gradient descent method. At first, newly added weak classifiers are trained based on the negative gradient information of the current model loss function, and then combines the trained weak classifiers into the existing model in an additive form. In short, the prediction model produced by GBR is an integration of weak prediction models. Support vector regression algorithm. Support vector regression (SVR) algorithm constructs hyperplanes or hyperplane sets in high or infinite-dimensional spaces, which can be used for classification regression or other tasks 36 . Intuitively, the farther the classification boundary is from the nearest training data point, the better the model because it can reduce the generalization error of the classifier. The SVR is trying to find a hyperplane that minimizes the distance from all data to the hyperplane.
Random forest regression algorithm. The random forest regression (RFR) is a classifier that contains multiple decision trees, and its output category is determined by the mode of the sum of all the individual trees 37 . RFR randomly uses variables and data of a data set to generate many classification trees that form a forest, and the prediction result is obtained by weighting the individual classification tree.

Results and discussion
Construction of data set. As shown in Fig. 6a, the dissipated energy of the UD-FRCP for various fiber contents are summarized based on the 185 impact simulation cases. It can be seen that the impact energy dissipation of the UD-FRCP increases from 2.83 to 7.14 J with increasing the fiber content from 20.34 to 39.92%, which indicates that the properties of UD-FRCP are rising with the more fiber content. The dissipated energy for the fiber content of 42% is almost three times larger than that for the fiber content of about 24%. In addition, although the total dissipated energy of the UD-FRCP is dispersed, we can get the tendency from the result. The dispersion of normalized dissipated energy, obtained by using FEA results divide the corresponding interval average value, is diminishing with increasing the fiber content as shown in Fig. 6b, where the normalized dissipated energy for the fiber contents in the range of 18 ~ 22%, 28 ~ 32%, and 38 ~ 42% are given. For low fiber contents of 18 ~ 22%, the data divergence can be as high as 45%, while it decreases to 20% for relatively high fiber contents of 38 ~ 42%. The divergence of dissipated impact energy of the UD-FRCP for the same fiber content should be ascribed to the difference of two-point function values of microstructures introduced by the uncertainty of fiber bundle positions in the cross-section. For low fiber contents, high uncertainty of fiber bundle positions leads to a noticeable difference of two-point correlation function values and finally results in an increased divergence of dissipated impact energy. Therefore, besides the fiber content, it is vital to consider the effects of microstructures on the impact protective performance of composites.
The sensitivity of parameters of the training algorithms. The model is trained by the GBR, SVR, and RFR algorithms in Scikit-learn 38 to study the difference of training algorithms and determine the appropriate algorithm for the present problem. We randomly choose 175 cases as a training data set during training the model and the rest 10 cases as the testing set. According to the definition of the No Free Lunch Theorem, it is known that no algorithm and parameter is suitable for all data sets 39 . Therefore, it is crucial to obtain the opti- www.nature.com/scientificreports/ mized parameters for a given algorithm to predict better. Consequently, the sensitivity analysis of the critical parameters for GBR, SVR, and RFR algorithms are performed. For the GBR algorithm, the fundamental parameters include the learning rate, the number of estimators, and the maximum depth of a single decision tree. In the present study, the learning rate is determined as 0.001 to ensure that the optimal parameters can be obtained while minimizing the computation cost. We mainly focus on the effects of the number of estimators and the maximum depth of a single decision tree on the training results of the GBR algorithm. For the RFR and the SVR algorithms, the sensitivities of the number of estimators and the kernel functions are investigated, respectively.
The sensitivities of the number of estimators and the maximum depth of a single decision tree on the training results of the GBR algorithm are shown in Fig. 7a,b, respectively. From Fig. 7a, it can be seen that the average error decreases quickly from 15.47 to 6.94%, and the maximum error decreases abruptly from 51.16 to 12.69% with increasing the number of estimators from 10 to 300. After that, the average and the maximum errors approach the plateaus at 6.94 and 12.69% as the number of estimators continues to increase. Generally, the number of estimators is equivalent to how many weak classifiers are used in the model. When the number of classifiers reaches a critical value, the accuracy of the training model will be saturated, resulting in the plateaus of errors with the continual increase of the number of estimators, as observed in Fig. 7a. Therefore, the number of estimators is determined as 300 for the GBR algorithm. As shown in Fig. 7b, the effect of the maximum depth of a single decision is non-monotonic. The average error decreases from 9.76 to 6.94%, and the maximum error decreases rapidly from 36.82 to 12.69% with increasing the maximum depth of a single decision tree from 2 to 3, then both the average and the maximum errors increase slightly with continually increasing the maximum depth of a single decision tree to 4. It can be well understood by the fitting ability of the GBR with increasing the maximum depth of a single decision tree. When the depth of a single decision tree is too small, the characteristics of the data set cannot be fitted well, whereas if it is too large, excessive parameters are introduced in the GBR algorithm, leading to unsatisfied fitting results due to the limited training data set, which we call over-fit 40 . As a result, the appropriate depth of the single decision tree is determined as 3 for the GBR algorithm in the present study.
The effect of the number of estimators on the training results of the RFR algorithm is also depicted in Fig. 7a. Similar to the GBR algorithm, the average and the maximum errors for the RFR algorithm decrease with increasing the number of estimators, and then they gradually decay to stable values at 7.87 and 17.26%, respectively, for the number of estimators more than 300. Therefore, the number of estimators is also determined as 300 for the RFR algorithm. The influences of various kernel functions that are radial basis function (RBF), polynomial function (Poly), and sigmoid functions (Sigmoid), on the training results of the SVR algorithm, are shown in Fig. 7c  www.nature.com/scientificreports/ 98.32% for Sigmoid, respectively. It can be seen that the training results with Poly are superior to those with the RBF and Sigmoid for the SVR algorithm, and therefore the Poly is used for the prediction by the SVR algorithm. The prediction results from the GBR, the RFR, and the SVR algorithms are given in Fig. 7d, and the average and the maximum errors are listed in Table 2. It is clearly shown that the GBR algorithm is more appropriate for the present study in terms of the average and the maximum errors when compared to the RFR and the SVR algorithms. The GBR algorithm gives the best predictions with an average error of 6.94% and a maximum error of 12.69%, as listed in Table 3, followed by the RFR algorithms with an average error of 7.87% and a maximum error of 19.61%. The SVR algorithm arises the worst predictions with an average error of 18.49% and a maximum error of 58.77%, which is several times larger than the results of the GBR algorithm and unacceptable for the present study. Such a significant error of the SVR algorithm should be ascribed to the intrinsic incapability of support vector algorithms as the SVR algorithm for the non-continuous regression problems, as in the present study. For a support vector algorithm, a data set is classified on a multi-dimensional plane to find the hyperplane and the maximum distance. If it is applied to a non-continuous regression problem, it is difficult to find the hyperplane, and the results of each category are voted to finally obtain an average value, leading to a relatively large error. In contrast, the GBR and the RFR algorithms in the Scikit-learn package are based on decision trees, which are suitable for dealing with non-continuous regression problems as employed in the present study. Specifically, for the RFR algorithm, the data set is randomly sampled into the decision trees, and the prediction result is averaged by all the votes of the multiple decision trees. While for the GBR algorithm, a first randomly selected decision tree is used to fit the data set, and the training error is then substituted into the second chosen decision tree randomly, and so on. With this method, the training result finally approaches the actual value step by step. As a result, both the GBR and the RFR algorithms can give good predictions, and the GBR algorithm is more appropriate for predicting the impact of the proactive performance of UD-FRCP. Further ML study will be performed to predict the dynamic mechanical behavior of other types of composites, such as two-or threedimensional fabric composites, in the future.
The effect of the number of cases of the GBR algorithm. In addition, we explored the effect of the number of training cases on the prediction error. We chose the training cases and predicting cases randomly, and ensured the number ratio of predicting cases to training cases almost constant for comparison. It can be seen that the average and maximum errors basically decrease as the training cases increases as given in Table 4, and an acceptable accuracy is reached for 175 training cases and 10 predicting cases.
The discussion on the uncertainty of the prediction. The uncertainties of the data set can be categorized into two aspects: aleatoric uncertainty and epistemic uncertainty when constructing the prediction model by using a finite size training set 41 . Aleatoric uncertainty is determined by the quality of the data set and the essence of the prediction task. There are inevitable measurement errors resulted from the precision of instruments while collecting data. Epistemic uncertainty is coming from the training process for the reason that it has intrinsic error between the prediction value and the true value due to the limitation of data size. The error will decrease with increasing the amount of training data.
In this work, the relationship between the microstructure and the penetration resistance of UD-FRCP was investigated. The data set was built by extracting the characteristics of the microstructure as input data. The absorbed energy of the composite after perforation were used in the FEA as output data. The error caused by anthropogenic factors is non-existent because the same simulation conditions are used for all of the cases.  www.nature.com/scientificreports/ Therefore, it is relatively reliable in terms of aleatoric uncertainty. For the epistemic uncertainty in this work, with the increase of training cases, the training data with random 175 cases are close to the distribution of the task, and the machine learning model can predict a value close to the true value with a maximum error of 12.69%. Therefore, it is also acceptable in terms of epistemic uncertainty. In summary, the machine learning model in this paper can effectively predict the penetration resistance of UD-FRCP by microstructure recognition.

Conclusions
The ML method is employed to predict the impact protective performance of UD-FRCP for the first time. The main conclusions can be summarized as follows.
1. A ML model is developed based on the 185 micro-scale simulation results to predict the macroscopic impact protective performance of UD-FRCP. The average error is 6.94%, and the maximum error is 12.69%, which is acceptable for ballistic impact problems of composites, providing a high-efficient method for analyzing the dynamic behavior of UD-FRCP. 2. The sensitivities of the critical parameters of training algorithms are investigated. The results showed that the increase in the number of estimators for GBR and RFR could increase the accuracy of the model, and the critical number of 300 is observed for the saturation of the accuracy. The increase of the maximum depth of a single decision tree for the GBR algorithm shows non-linear effects on the accuracy of the model. 3. The decision trees based algorithms are appropriate for the present study problem, and the GBR algorithm is found to give the best prediction results. More ML models will be built for other types of composites in the future. 4. The number of 185 cases is appropriate for this problem compared to other numbers of cases. The accuracy is guaranteed without a very large computational effort.