Interpretable Machine Learning Models for Punching Shear Strength Estimation of FRP Reinforced Concrete Slabs

: Fiber reinforced polymer (FRP) serves as a prospective alternative to reinforcement in concrete slabs. However, similarly to traditional reinforced concrete slabs, FRP reinforced concrete slabs are susceptible to punching shear failure. Accounts of the insufficient consideration of impact factors, existing empirical models and design provisions for punching strength of FRP reinforced concrete slabs have some problems such as high bias and variance. This study established machine learning ‐ based models to accurately predict the punching shear strength of FRP reinforced concrete slabs. A database of 121 groups of experimental results of FRP reinforced concrete slabs are collected from a literature review. Several machine learning algorithms, such as artificial neural network, support vector machine, decision tree, and adaptive boosting, are selected to build models and com ‐ pare the performance between them. To demonstrate the predicted accuracy of machine learning, this paper also introduces 6 empirical models and design codes for comparative analysis. The com ‐ parative results demonstrate that adaptive boosting has the highest predicted precision, in which the root mean squared error, mean absolute error and coefficient of determination of which are 29.83, 23.00 and 0.99, respectively. GB 50010 ‐ 2010 (2015) has the best predicted performance among these empirical models and design codes, and ACI 318 ‐ 19 has the similar result. In addition, among these empirical models, the model proposed by El ‐ Ghandour et al. (1999) has the highest predicted accuracy. According to the results obtained above, SHapley Additive exPlanation (SHAP) is adopted to illustrate the predicted process of AdaBoost. SHAP not only provides global and indi ‐ vidual interpretations, but also carries out feature dependency analysis for each input variable. The interpretation results of the model reflect the importance and contribution of the factors that influ ‐ ence the punching shear strength in the machine learning model.


Introduction
Reinforced concrete slabs are one of the common horizontal load-carrying members in civil engineering, and widely applied in bridges, ports and hydro-structures. Since there is no beam in the flat slab under longitudinal load, the punching failure of reinforced concrete slab occurred easily. Many researchers have carried out numerous experiments on the punching shear resistance of reinforced concrete slabs, and obtained successful results [1][2][3]. However, steel bars are prone to corrosion, which will result in the shortening of the actual service life [4]. In recent years, with the development of technology, the durability of structure has attracted people's attention. For coastal areas and the areas of using chlorides such as deicing salt, the actual service life of structures is often much lower than their design service life, resulting in massive losses [5].
Fiber reinforced polymer (FRP) is a material which has many advantages such as light, high strength and corrosion resistance. In a corrosion environment, to solve the problem of short actual service life of structure, FRP bar can be applied as an alternative to steel bars in concrete structures. This is because FRP bar can be appropriate for service demands of concrete structure in severe environments which contain plenty of chloride ion and sulfate ion. However, FRP bar has a low elastic modulus, which will result in FRP reinforced concrete slabs being more prone to punching failure than reinforced slabs (Figure 1). Many experimental studies show that [6][7][8][9][10][11] under the same conditions, the punching shear strength and stiffness of the cracked FRP reinforced concrete slabs decrease faster than reinforced concrete slabs. Matthys and Taerwe [6] concluded that the crack development of slabs and brittleness of punching failure were significantly affected by the bond performance of FRP grid reinforcement through experimental results. The experimental results of Ospina et al. [7] indicated that the bond behavior between FRP bar and concrete have a significant impact on punching shear capacity of FRP reinforced concrete slabs. It is shown that the slab thickness and the concrete compressive strength were of considerable influence on the punching shear capacities of FRP reinforced concrete slabs by Bouguerra et al. [9]. Ramzy et al. [12] reported that the punching shear behavior of flat slab was related to the size effect of construction and the Young's modulus of FRP reinforcement. In terms of theoretical models, most of the punching shear strength computational formulas of FRP reinforced concrete slabs were derived from traditional reinforced concrete flat with modifications to account for FRP [13]. Some current design specifications such as the GB 50010-2010 [14] and the ACI 318-14 [15] regard the eccentric shear stress model as theoretical basis. Based on the ACI 318-11, El-Ghandour et al. [16] suggested that one considers the impact of elastic modulus of FRP bars, and come up with an improved equation for punching shear strength. After that, El-Ghandour et al. [17] used the same method to modify the design formula of the BS 8110-97. Matthys and Taerwe [6] considered the effect of the equivalent reinforcement ratio, and proposed a modification of the BS 8110-97. On this basis, Ospina et al. [7] proposed an empirical model for computing the punching shear strength of FRP slabs by modifying the relation between the punching shear capacity and the equivalent reinforcement ratio. On the basis of the probability of exceedance, Ju et al. [18] proposed a new approach for analyzing the punching shear strength of FRP reinforced two-way concrete slabs by using Monte Carlo simulations. However, the aforementioned empirical models adopted some simplifications during theoretical derivations, thus the empirical models were unable to consider all of the influential factors. What's more, the parameters in the aforementioned empirical models were determined by the traditional regression analyses from experimental results. Therefore, the accuracy of the models is highly dependent on the choices of theoretical models and quality of the databases.
In recent years, with the development of artificial intelligence, some algorithms with data at the core have emerged [19]. Among these algorithms, machine learning has received remarkable attention of researchers, and there have been many successful examples [20][21][22][23][24]. In structure engineering, Hoang et al. [25] constructed machine learning based alternatives for estimating the punching shear capacity of steel fiber reinforced concrete (SFRC) flat slabs. Hoang et al. [26] presented the development of an ensemble ma-chine learning model to predict the punching shear resistance of R/C interior slabs. Mangalathu et al. [27] build an explainable machine learning model to predict the punching shear strength of flat slabs without transverse reinforcement. In addition, some researchers [28,29] even use the atomistic simulations as the input parameters of machine learning to predict the performance of materials and structures, which has also seen success. To the best of the authors' knowledge, however, no study examined the interpretable machine learning models in predicting the punching shear strength of FRP reinforced concrete slabs. Therefore, this study intends to fill this research gap. The deep relationship between material properties and punching shear strength will be found if the machine learning research on FRP reinforced concrete slabs is carried out. It is helpful to refine the traditional empirical models and design codes and make them more accurate than before. In this paper, an experimental database for the punching shear strength of FRP reinforced concrete slabs is first compiled, and used for training, validating, and testing machine learning models.
Four machine learning algorithms, namely artificial neural network (ANN), support vector machine (SVM), decision tree (DT) and adaptive boosting (AdaBoost), are employed for punching shear strength prediction. Then, by comparing the performance of machine learning models, their efficiency and accuracy can be determined. These comparisons are valuable for identifying the efficiency and prediction ability of the machine learning models. However, there is problems which need to be solved; for example, the values of hyper-parameters in machine learning are often difficult to determine. Therefore, to assist the improvement of the predicted performance of machine learning models, optimization methods such as particle swarm optimization (PSO) and empirical method will be used. In addition, since previous models based on machine learning found it difficult to explain the predicted mechanism, the predicted results of models were somehow unconvincing [27]. This paper uses SHapley Additive exPlanation (SHAP) [30] to explain the predicted result of model. Distinct from other machine learning papers explained by SHAP, the kernel explainer will be used in this paper, which is appropriate for all the machine learning models. Not only can SHAP carry out analysis of feature importance for input factors, it can also determine whether the impact of input features on predictions is positive or negative. In addition, SHAP can make researchers realize the interrelation between input features, and how each input feature will influence the final predicted value for a single sample. Consequently, the emergence of SHAP renders the predicted results of machine learning more convincing than before.

Experimental Database of FRP Reinforced Concrete Slab
To ensure the accuracy of predicted results, a database with adequate samples is required for model training. In this paper, 121 groups of experimental results [6][7][8][9][10]12,[31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] of FRP reinforced concrete slabs under punching shear tests were collected, and made it into a database. In the data set, 80% of the data was used as a training set (10% of this was used as validation set), and remaining 20% were used as test set. In this study, 6 critical parameters are selected to characterize the punching shear strength of FRP slabs such as the types of column section, cross-section area of column (A/cm 2 ), slab's effective depth (d/mm), compressive strength of concrete (f'c/Mpa), Young's modulus of FRP reinforcement (Ef/Gpa), and reinforcement ratio (ρf/%). According to the findings of Ju et al. [18], these 6 parameters have a great influence on punching shear strength of FRP slabs. Therefore, these 6 parameters are considered as input parameters and the punching shear strength of FRP slabs is considered as an output parameter of machine learning algorithms. To simplify the names of the 6 input parameters and the 1 output parameter, we used x1 to x6 and y to represent them, respectively. The column section has 3 types: square (x1 = 1), circle (x1 = 2), rectangle (x1 = 3). The distribution of the dataset is shown as Table 1 and Figure 2, and the detailed information is shown as Appendix A.  Before starting the machine learning procedure, the data of input variables should be preprocessed. In this paper, we used deviation standardization as the preprocessing method; this can be written as: where x is the sample before feature scaling; x * is the sample after feature scaling; m is the total number of samples. Since the sample x contains features of several input parameters, it can also be named the feature vector.

Machine Learning Algorithms
Generally, the implementation of machine learning has four stages: (a) divide the database into training set and test set; (b) apply the training set to model training; (c) check whether the accuracy requirements are met; d) output the predicted model for test or adjust the values of hyper-parameters. The flowchart for this procedure is shown as Figure  3. A total of 4 machine learning models, namely, ANN, SVM, DT and Adaboost, are selected in this study. As for ANN, an input layer, hidden layers and an output layer are formed for the predicted results, where several neurons are applied. As for SVM, a prediction equation is formed based on the predicted output in which the error should be smaller than a fixed value. As for DT, the whole database is separated into tree-like decisions which are based on one or several input characteristics, and the output is the tree which falls into the data. As for the AdaBoost, a group of weak learners are collected to form a strong learner which achieves better predicted results. Notably, the predicted performance of model has relations with the value of hyper-parameters. Therefore, in this paper, some methods will be selected to better determine the value of hyper-parameters.

Artificial Neural Network
Artificial neural network (ANN) is a supervised learning algorithm that was originally used to simulate the neural network of human brain [50]. Generally, a simple ANN model consists of an input layer, one or two hidden layers, and an output layer. There are several neurons in each layer of ANN, where the number of neurons in the input layer and output layer depends on the number of input parameters and output parameters, respectively. In this paper, the number of neurons in the input layer and output layer can be seen in Table 1, which equals 6 and 1, respectively.
Except for the output layer, the neurons of each layer needs to be activated by activation function. In ANN, generally, the sigmoid function is used as activation function [51], which expression can be written as: where b is the bias unit; w denotes the weight vector of sample. According to Rumelhart et al. [52] described, when ANN is used for regression forecasting, the cost function of model as follows: where ypred denotes the predicted value of machine learning; C is the regularization coefficient.
Since the cost function represents the error between the actual value and the predicted value, to improve the accuracy of predicted result, a method similar to gradient descent, in order to decrease the value of cost function, is required. The gradient descent is expressed as: where α is learning rate. Before starting the training procedure, there are some values of hyper-parameters need to be set. For neuron number of hidden layer, Shen et al. [53] mentioned a method to determine their value range, which method as follows: where X is the number of neurons in the input layer; H is the number of neurons in the hidden layer; O is the number of neurons in the output layer; a is the constant between 1 to 10. According to the Equation (6), it can be ascertained that the range of neurons in hidden layer ranges from 7 to 13. In order to better compare the performance of models which have different neuron numbers in the hidden layer, 3 indicators are introduced, which names root mean squared error (RMSE), mean absolute error (MAE) and coefficient of determination (R 2 ), which are, respectively, defined as: The performance of ANN which have different neuron numbers in hidden layer is shown as Table 2. According to the comparative results, the number of neurons in the hidden layer H is determined as 8. In addition, the values of other hyper-parameters confirmed by manual adjustment are: the number of hidden layers equal 1; the training times equal 5000; the learning rate α equals 5 × 10 −5 ; the regularization coefficient C equals 20.

Support Vector Machine
Support vector machine (SVM) was originated from statistical learning theory, and it was firstly proposed by Cortes and Vapnik [54] in 1995. When SVM is used for regression analysis, it can also be called support vector regression (SVR). The predicted method of SVR maps the data to high-dimensional space and use hyperplane for fitting, and these data are generally hard to fit into low-dimensional space [55]. The expression of SVR reflect the relation between input factors and predicted value [56], it can be written as: where ϕ(x) means feature vector after mapping. With the same as ANN, SVR also has the cost function to represents the performance of model. The cost function can be written as: where ξ and ξ * are the slack variables, they represent the errors of the sample exceeding the lower and upper bounds of the hyperplane, respectively.
To minimize the value of cost function, Lagrange multiplier method is used for transforming the original problem to its dual problem. The Lagrange multiplier method is formulated in the following form: where μ, μ * , λ, λ * denote the Lagrange multiplier; ε is the maximum allowable error. After transforming, the expression of dual problem as follows: The Karush-Kuhn-Tucker (KKT) conditions need to be met for Equation (13), which can be written as: According to Equations (13) and (14), the weight vector w and bias unit b can be acquired. Finally, the expression of SVR as follows: where κ(xi, xj) denotes kernel function, which can used to express the inner product of eigenvectors in high-dimensional space. The kernel function has many types, and the radial basis function (RBF) kernel function is selected in this paper, which can be expressed as: where γ > 0 is the coefficient of RBF kernel function, which is a hyper-parameter. Except that, in this paper, the values of other hyper-parameters in SVR need to be determined: the maximum allowable error ε, the regularization coefficient C. There are many methods which can be used for numerical adjustment of hyper-parameters, in this paper we have selected particle swarm optimization. Particle swarm optimization (PSO) is an optimized algorithm inspired by the foraging behavior of bird flock [57,58]. On account of the fact that PSO has a large number of advantages, such as being easily comprehendible and containing, it is usually used for hyper-parameter adjustment in SVM [59,60]. For every particle in PSO, there is two parameters to determine the search result, which are position p and velocity v, respectively. These 2 parameters are defined as [58]: where t is the number of iterations; d * denotes the dimension of search space; W is the inertial coefficient; C1 and C2 are the self-learning factor and the global learning factor; r1 and r2 are random number between (0, 1); P is the location with the best fitness among all the visited locations of each particle; G is the location with the best fitness among all the visited locations of all the particles. In this paper, some hyper-parameters in Equations (18) and (19) have been confirmed: the maximum of iteration equals 150; the dimension of search space d * equals 3; the inertial coefficient W equals 0.8; the total number of particles equals 40; the self-learning factor C1 and global learning factor C2 are all equal to 2. The hybrid machine learning approach consists of PSO and SVR can be called PSO-SVR, whose best value of hyperparameters also can be determined by calculating: the maximum allowable error ε equals 10; the regularization coefficient C equals 4238.58; the coefficient of RBF kernel function γ equals 2.88.

Decision Tree
Decision tree (DT) served as a supervised learning algorithm that is proposed by Quinlan [61]. In the beginning, DT only could be used for conducting the classification forecasting, and its types less such as iterative dichotomiser 3 (ID3) and classifier 4.5 (C4.5). In this paper, classification and regression tree (CART) [62] is selected to make a study of regression prediction. The dividing evidence of CART is the variance of samples, which can be written as: where R1 and R2 denote the sample set after dividing; c1 and c2 are mean of samples in R1 and R2, respectively. With the division of sample set, gradually, the purity of leaf nodes will improve.
Pruning is a method to avoid the over-fitting in DT, and it contains two basic methods: pre-pruning and post-pruning [63,64]. In this paper, the specific way of pruning is to adjust the maximum of depth, which equals 5.

Adaptive Boosting
Adaptive boosting (AdaBoost) is one of the ensemble learning algorithms that can combine multiple weak learners into strong learner [65]. To be distinct from aforementioned machine learning algorithms, each sample in AdaBoost has a subsampling weight which will constantly adjust [66]. At the beginning of the training process, the weights of samples need to be initialized [67], which can be written as: where w * ki denotes the weight of the i-th sample in the k-th weak learner. The largest error of the k-th weak learner on the training set as follows: where Gk(xi) denotes the predicted value of the i-th sample in the k-th weak learner. The relative error eki of each sample and the error ratio ek of the k-th weak learner can be written as: After that, the weight βk of the k-th weak learner can be calculated, the value of which will increase as the error ratio ek decreases. The weight βk can be defined as: Afterwards, the error ratio of the next weak learner can be reckoned, and the weights of samples need to be updated: where Zk is the normalization factor. Finally, the expression of the strong learner can be written as: where K is the total number of weak learners; g(x) is the median of all the βkGk(x). The algorithm of the weak learner needs to be confirmed, and PSO-SVR is chosen in this paper.

Predicted Results
The result comparison between ANN, PSO-SVR, DT and AdaBoost are provided in Table 3. The four machine learning models shows high predicted accuracy reflected in R 2 . Especially AdaBoost and PSO-SVR, the predicted accuracy of the former is highest among these four AI models in training set, and the latter is in test set. It is noted that although the forecasting performance of ANN for the test set is better than DT and just slightly lower than PSO-SVR and AdaBoost, the performance for training set is much lower than the other three models.  Figure 4, it can be found that the fitting condition of ANN is not very suitable when the true value is larger than 1400 kN. Maybe this situation is related to insufficient data of punching shear strength larger than 1400 kN. Therefore, it is noted that if we want to use this ANN model for predicting the punching shear strength of FRP reinforced concrete slabs, the value range of the predicted object should be restrained under 1400 kN. In addition, from the figure, it can be seen that ANN focuses more on fitting samples which are located in the interval with dense sample distribution. Consequently, a way to improve the overall predicted accuracy of ANN is to ensure that the samples are evenly distributed. It is demonstrated in Figure 5 that PSO-SVR shows the good predicted performance. Although there some predicted values have a relatively larger error between the actual values, the overall predicted accuracy is higher than ANN and DT. Figure 6 depicts the predictions of DT on the dataset; there doubt that the predicted result of DT for the test set is worse than ANN and PSO-SVR. In training set, the predictions of DT are almost no error with the punching shear strength is larger than 1200 kN, it is related to the predicted mechanism of DT [68]. Moreover, the situation also can be seen from the predicted results of DT on punching shear strength less than 1200 kN, the predicted value of which is same when the punching strength is in the identical range. Meanwhile, in test set, DT has a low predicted accuracy for samples with punching shear strength of around 1200 kN, i.e., the DT model proposed in this paper has the risk of over-fitting when the punching shear strength larger than 1200 kN. Obviously, no matter ANN or DT, the predicted accuracy of them is closely related to the number of samples in each value range. Furthermore, in this paper, it can be concluded that neither ANN nor DT has mined the true intrinsic connection between the data. Figure 7 shows the predicted results of AdaBoost for the training set and the test set, it can be seen that all of the dots are closely distributed around the best fitting line. Since PSO-SVR is chosen as the weak learner of Ada-Boost, some wildly inaccurate predicted values are improved, and some comparatively accurate predicted values are retained. Moreover, compare Figure 4 to Figure 7, it is no hard to see that the fitting situation of AdaBoost reveals the best predicting performance.

Comparison with Traditional Empirical Models
In order to reflect the accuracy of the machine learning calculation results, this paper introduces design specifications and empirical models proposed by researchers. In this study, 3 punching shear strength models are selected from design codes and reviewed: GB 50010-2010 (2015) [14], ACI 318-19 [69] and BS 8110-97 [70]. Most design codes adopt the same function for both FRP slabs and reinforced concrete slabs. The punching shear strength of slabs as a function of the concrete compressive strength. Some design codes also take the reinforcement ratio, size effect, depth of the slab and so on into consideration. Three design specifications and three empirical models are adopted in comparative analysis between calculations. The formulas of selected traditional empirical model are shown in Table 4.
βh is the sectional depth influence coefficient; ft is the design value of tensile strength; b0,0.5d is the the perimeter of the critical section for slabs and footings at a distance of d/2 away from the column face; βs is the ratio of the long side to the short side of the sectional shape under local load or concentrated reaction force; αs is the influence coefficient of column type; fcu is the compressive strength of concrete cube; b0,1.5d is the perimeter of the critical section for slabs and footings at a distance of 1.5d/2 away from the loaded area; Es is the Young's modulus of steel bar. Table 5 summaries experimental results of machine learning models and traditional models. Obviously, the predicted performance of AdaBoost is better than other predicted models reflected in RMSE, MAE and R 2 . ANN, DT and PSO-SVR, the other three AI models, although their predicted results are not best, they show a better performance compared with traditional empirical models. Among these traditional models, formula 1 (GB 50010-2010 (2015)) has the highest predicted accuracy, the RMSE, MAE and R 2 of which are 150.41, 98.48 and 0.73, respectively. Notably, the predictions between formula 1 and formula 2 only has a little error; perhaps it is related to the fact that they are both derived from the eccentric shear stress model [14,69]. In addition, the forecasting performance of the model proposed by El-Ghandour et al. is better than the other traditional models proposed by other researchers. In brief, the RMSE and MAE of traditional empirical models are generally greater than 100, and the R 2 of them ranged from 0.6 to 0.8. However, the RMSE and MAE of machine learning models are less than 100, and the R 2 of them are generally greater than 0.9 and even close to 1.0. Consequently, it can be said that the difference in the predicted accuracy between traditional empirical models and machine learning models is obvious.  Figure 8 shows the predicted results of traditional empirical models and machine learning models for the whole data set. To facilitate comparison between different types of models in these figures, the red dots are used to represent the predicted result of national codes. Afterwards, the blue and purple dots are used to represent the predicted results of traditional models proposed by researchers and machine learning models, respectively. It can be found that the predicted error of traditional models is relatively large, and the machine learning models reach the opposite conclusion. The forecasting result of Figure 8a is similar to Figure 8b, and this conclusion is also reflected in Table 5. In addition, although the predicted accuracy of formula 1 and formula 2 is better than other traditional empirical models, their predicted results are still unsafe. On the contrary, the predicted results of formula 3 to formula 6 tend to be conservative even though they do not have the best predicted performance. From the predicted performance of the machine learning models for the complete database, it can be found that the predicted results of ANN and PSO-SVR are likely to be unsafe when the punching shear strength smaller than 200 kN. On the side, when the punching shear strength larger than 1200 kN, the forecasting result of ANN is conservative. Furthermore, except formula 4, the predicted accuracy of all the empirical models are better than ANN when the punching strength larger than 1400 kN. Consequently, it can be concluded that ANN is very dependent on samples.

Shapley Additive Explanation
In this study, SHapley Additive exPlanation (SHAP) is applied to explain the predicting results given by the AdaBoost model. SHAP originates from game theory and is an additive model, i.e., the output of the model is a linear addition of the input variables [71]. The expression of SHAP can be defined as: where ybase is the baseline value of model; n is the total number of input variables; f(xij) is the SHAP value of the xij. According to Equation (29), it can be concluded that the predicted value of each sample depends on the SHAP value of each feature in the sample. In other words, each feature is a contributor for final predicted value. The SHAP value can be positive or negative, depending on the influence trend of each feature to predicted result. SHAP contains several explainers, in current machine learning papers it often uses the tree explainer to illustrate its models, which is only suitable for a part of machine learning algorithms. Consequently, the kernel explainer will be used to illustrate the machine learning model which has the best predicted performance in this paper.

Global Interpretations
According to the experimental results obtained above, the predicted results of Ada-Boost, which has the best predicted performance, was explained. The predicted results of AdaBoost can be interpreted by SHAP in different ways. Firstly, the feature importance analysis of input parameters is shown in Figure 9a. It can be seen that the importance of each input parameter is non-negligible, that is each input parameter will have a different degree of influence on the predicted results. This is consistent with the existing experimental conclusions [4,72]. According to the analysis results, it can be understood that the SHAP value of each parameter are x1 equals 53.45, x2 equals 44.50, x3 equals 126.36, x4 equals 21.57, x5 equals 17.23, x6 equals 52.96, respectively. Obviously, the slab's effective depth (x3) has the most critical influence on the punching shear capacity, more than twice as much as types of column section (x1). Moreover, though the Young's modulus of FRP reinforcement (x5) has the least importance on the predicted results, sometimes ignoring this feature will cause unnecessary errors in the predicted results of model. Consequently, using these 6 features as the input parameters of machine learning algorithm is reasonable.  Figure 9b is a SHAP summary plot of the features, which displays the distribution of the SHAP values of each input variable in each sample and demonstrates the overall influence trends. In the figure, the x-axis is the SHAP value of input parameters, the y-axis is the input parameters, ranked by importance. The points in the figure represent the samples in the dataset and their colors represent the feature value, for which color from blue to red represents the value from small to large. From the figure, it can be seen that a few of SHAP values of samples in x3 are very high, which are around 600. It can be said that in these samples, the value of x3 has decisive effect on the predicted results. In these input variables, the slab's effective depth (x3), types of column section (x1), reinforcement ratio (x6), cross-section area of column (x2), Young's modulus of FRP reinforcement (x5) all have the positive influence on the punching shear strength. However, the influence trend of the compressive strength of concrete (x4) on the predicted results is hard to definite, which is very complex.

Individual Interpretations
Except the global interpretations, SHAP also can provides the individual interpretation for each sample in the dataset. In this study, 2 samples are selected, the predicted process of which are illustrated as Figure 10. In the figure, the gray line represents the baseline value, which equals 403.62, and the colorful lines represents the decisive process of each feature. In Figure 10, it can be seen that from the baseline value, every feature will change the value to a different degree, and finally predicted values 419.27 and 501.00, respectively, will be obtained. In the predicted process of the 2 samples, types of column section (x1), slab's effective depth (x3), the compressive strength of concrete (x4) have negative influence on predicted results; whereas, the cross-section area of column (x2), Young's modulus of FRP reinforcement (x5), reinforcement ratio (x6) have positive influence on it. Furthermore, reinforcement ratio (x6) acts as determinant in forecasting process of sample 1, and the cross-section area of column (x2) plays an essential role in sample 2. In summary, both the global interpretations and the individual interpretations give a detailed insight into the process by which each input feature affects the final predicted value.

Feature Dependency
In addition to explaining the predicted mechanism of model, SHAP can also reveal the interaction between the specific feature and its most closely related feature, which is shown in Figure 11. It can be seen from the Figure that with the values increasing of x1, x2, x3, x6, their SHAP values will also be increasing, but the relationships between x4, x5 and their SHAP values are complicated and nonlinear. As the values of x4, x5 increase, the SHAP values of them will increase or decrease; thus, it is difficult to describe the concrete relationships. However, the one that interacts most closely with x1, x2, x5, x6 is x3, the features most closely related to x3 and x4 are x2 and x6, respectively. According to the results of feature dependency, further experiments can be carried out to analyze the relationship between the two variables that interact most closely, and analyze the collaborative impact of these two variables to output. It might be helpful to improve the punching shear resistance of FRP reinforced concrete slab-column connections, and the traditional empirical models cannot accomplish this because the traditional models, relatively, have the worse predicted performance.

Conclusions
In this paper, several machine learning algorithms, ANN, SVR, DT and AdaBoost, are adopted for punching shear strength of FRP reinforced concrete slabs. A total of 121 groups of experimental tests are collected and the test variables are divided into 6 inputs (types of column section, cross-section area of column, slab's effective depth, compressive strength of concrete, Young's modulus of FRP reinforcement, reinforcement ratio) and 1 output (punching shear strength). Random subsets of 80% and 20% of the data were used for training and testing, respectively, and 10% of the data in training set were used for validating. In addition, 3 indicators (RMSE, MAE, R 2 ) are introduced to assist the predicted performance. For improving the predicted performance, the empirical method and PSO are utilized for adjusting the values of hyper-parameters in ANN and SVR, respectively. After tuning, both ANN and SVR display the good predicted performance. According to the experimental results, the best performance model is AdaBoost, the RMSE, MAE, R 2 of which are 29.83, 23.00, 0.99, respectively. The predicted performance of ANN, DT and PSO-SVR are worse than AdaBoost, but they are all better than the 6 traditional empirical models introduced in this paper. The predicted result of ANN for training set has a relatively large deviation when the punching shear strength exceeding 1400 kN, but DT has almost no error. Among these traditional models, GB 50010-2010 (2015) has the least deviation reflected in RMSE and R 2 , the values of them are 150.41 and 0.73, respectively. Moreover, the model proposed by El-Ghandour et al. (1999) has the best forecasting performance among three traditional models proposed by researchers.
SHAP is used to interpret the predicted process of AdaBoost with the best performance in this study. On the basis of the result of global interpretations, it can be known that every input variable has the different degrees of influence on predicted results, which demonstrates the rationality of input variables selection. Furthermore, the slab's effective depth has the highest impact factor on the predicted results, equaling 126.36, more than twice as much as types of column section, which equals 53.45. In addition, the Young's modulus of FRP reinforcement has the lowest impact factor, equaling 17.23. From the SHAP summary plot, it can be understood that all input variables have positive influence on predicted results except the compressive strength of concrete, it is because its influence trend is relatively fuzzy. According to the individual interpretations, the forecasting process of single sample is shown, which can help researchers to understand the predicted process of model roundly. In addition, SHAP also provides the results of feature dependency, which is helpful to researchers in understanding the relationship between the specific feature and its most closely related feature. According to the analysis results of SHAP, the relationship between material properties and punching shear strength may really be revealed. Acknowledgments: This study is supported by the Engineering Research Centre of Precast Concrete of Zhejiang Province. The help of all members of the Engineering Research Centre is sincerely appreciated. We would also like to express our sincere appreciation to the anonymous referee for valuable suggestions and corrections.

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

Appendix A. Test Data of FRP Reinforced Concrete Slabs
The sources of the database and the specific material properties are listed in Table  A1. The input variables are the types of column section x1, cross-section area of column x2 (A/cm 2 ), slab's effective depth x3 (d/mm), compressive strength of concrete x4 (f'c/MPa), Young's modulus of FRP reinforcement x5 (Ef/GPa), and reinforcement ratio x6 (ρf/%), respectively. The output variable is the experimental value of punching shear strength y (V/kN). The column section has three types: square (x1 = 1), circle (x1 = 2), rectangle (x1 = 3). In order to use the experimental data for comparative analysis between models, the concrete compressive strength of columns with different cross-section types (cylinder and cube) was unified. Furthermore, to prevent errors caused by different data precision, retain two significant digits to the decimal point for input parameters except for column section type x1.