Comparative Analysis of Machine Learning Models for Predicting Interfacial Bond Strength of Fiber‐ Reinforced Polymer‐Concrete

1. Introduction

Fiber-Reinforced Polymers (FRPs) are essential in modern construction due to their high strength-to-weight ratio, corrosion resistance, and ease of application, making them ideal for retrofitting concrete and masonry structures. Retrofitting unreinforced masonry (URM) buildings in earthquake-prone regions with FRPs effectively improves seismic resistance, enhances load-bearing capacity, and increases ductility, addressing vulnerabilities caused by outdated construction techniques [1,2].

Traditional empirical models for predicting FRP-concrete bond strength have shown limitations in accurately capturing the complex interactions between these materials. This has led to a growing interest in applying machine learning (ML) techniques, which are capable of handling large datasets and uncovering intricate patterns that might be missed by conventional approaches. Recent studies have demonstrated the effectiveness of machine learning algorithms.

Wu and Jiang's research (2013) focused on quantifying and modeling interfacial bond parameters through both experimental and analytical studies. They developed a comprehensive database comprising 628 shear tests of externally bonded FRP joints. Their analysis identified key factors affecting bond parameters, particularly focusing on the bond strength and fracture energy. Additionally, they derived new models for the width factor as a function of width ratio and concrete strength, offering more accurate bond strength and fracture energy predictions. This work contributed to understanding how bond behavior is influenced by variables such as the width of the FRP and the strength of the concrete [3].

Zhou et al. (2020) developed an artificial neural network (ANN) model to predict the interfacial bond strength between FRP and concrete. Utilizing a substantial database of 969 test results, they trained the ANN model using backpropagation neural networks (BPNN), incorporating weighted values, biases, and transfer functions. Their model demonstrated higher accuracy than 20 existing models, with significantly lower predictive errors. The study also provided an explicit formula derived from the ANN model, which can be used in practical design applications for predicting bond strength [4].

In their research, Li et al. (2020) compiled a large dataset to assess the bond strength and bond stress-slip behavior of FRP-reinforced concrete members. Their study also incorporated an evaluation of environmental durability factors that could affect bond performance. The study is crucial for its contribution to understanding how reduction factors in bond strength are influenced by both the environment and material properties [5].

Su et al. (2021) explored the potential of machine learning approaches, including multiple linear regression, support vector machine (SVM), and artificial neural networks (ANN), for predicting the interfacial bond strength (IBS) between FRP and concrete. Two datasets were used: Dataset 1 (122 IBS values), and Dataset 2 (136 IBS values). The variables in Dataset 1 include FRP properties such as elastic modulus, tensile strength, and bond length, while Dataset 2 uses fewer input features but focuses on stiffness and groove dimensions. Among the machine learning models, the SVM model demonstrated superior predictive ability with R² values of 0.79 for Dataset 1 and 0.85 for Dataset 2. This research highlights the capability of advanced machine learning models in predicting FRP-concrete bond strength [6].

The study by Haddad and Haddad (2021) employed an artificial neural network (ANN) technique to predict the bond strength between FRP and concrete. The researchers compiled an extensive dataset comprising over 440 data points, allowing the ANN model to assess various parameters such as FRP properties, concrete compressive strength, and bond length. The model demonstrated impressive predictive accuracy, achieving high correlation coefficients (R² nearing 0.98), and significantly outperformed traditional empirical methods. This study underscores the effectiveness of ANN models in predicting complex interactions between FRP and concrete in structural applications [7].

Chen et al. (2021) developed a prediction model using the ensemble learning method known as Gradient Boosted Regression Trees (GBRT) to estimate the interfacial bond strength between FRP and concrete. Their model was trained using a comprehensive database of 520 tested samples. The GBRT model demonstrated superior accuracy, with R² values of 0.9627 during training and 0.9269 during testing. The model outperformed both ANN and SVM models, making it a robust tool for predicting FRP-concrete bond strength in practical engineering applications [8].

Barkhordari et al. (2022) introduced hybrid models that combine population-based algorithms (Bald Eagle Search, Manta Ray Foraging Optimization, and Runge-Kutta optimizer) with ANNs to estimate FRP-concrete interfacial bond strength. The study utilized a dataset of 969 experimental samples and found that the RUN-ANN model achieved the best performance, with an R² value of 0.92. Moreover, the Shapley Additive exPlanations (SHAP) method was used to interpret the model, identifying FRP bond length and width as key factors influencing bond strength predictions [9].

Alabdullh et al. (2022) applied a hybrid ensemble machine learning approach to predict the bond strength of FRP laminates bonded to concrete. Using a database of 136 samples, the team trained and validated several standalone machine learning models, including ANN, Extreme Learning Machine, and GPR. The hybrid ensemble model (HENS) achieved superior predictive accuracy, with an R² value of 0.9783 for training and 0.9287 for testing, demonstrating its effectiveness in overcoming overfitting issues common in traditional machine learning models [10].

In their study, Kim et al. (2022) employed ensemble machine learning techniques to predict FRP-concrete interfacial bond strength. Using a dataset of 855 single-lap shear tests, they evaluated various ensemble methods, with the CatBoost algorithm achieving the best performance. The algorithm showed a root mean square error of 2.310, a covariance of 21.8%, and an R² value of 96.1%, suggesting its high accuracy and potential for practical applications [11].

This research aims to advance existing methodologies by developing a robust machine learning model specifically designed to predict FRP-concrete bond strength, utilizing a comprehensive dataset of experimental results. The study evaluates several machine learning models to identify the most significant input variables and further interprets their predictions using advanced methods for feature importance analysis. Through this approach, the study aims to contribute to the ongoing efforts in optimizing FRP applications in civil engineering, ensuring more reliable and efficient reinforcement solutions.

2. Materials and Methods

2.1. Multiple Linear Regression

Linear regression analysis is a statistical method used to determine the causal effect of one or more independent variables on a dependent variable. When the problem involves one dependent variable ($Y$) and multiple independent variables ($X_1,X_2,\dots ,X_p$), multiple linear regression can be used. The model can be described by the equation (1):

$$Y={\mathsf{\beta }}_0+{\mathsf{\beta }}_1{X}_1+{\mathsf{\beta }}_2{X}_2+\dots +{\mathsf{\beta }}_p{X}_p+\mathsf{ϵ}$$

(1)

In this equation ${\mathsf{\beta }}_0$ is the intercept, ${\mathsf{\beta }}_1,{\mathsf{\beta }}_2,\dots ,{\mathsf{\beta }}_p$ are the coefficients of the independent variables, $ϵ$ represents the error term.

To estimate the parameters (${\beta }_0$, ${\beta }_1,{\beta }_2,\dots ,{\beta }_p$) of the linear model, it is used the least squares method. The goal is to minimize the sum of the squares of the differences between the observed values and the values predicted by the model. This is done using the following equation (2):

$$\widehat{\beta }={\left(X^{T}X\right)}^{-1}{X}^{T}Y$$

(2)

where $\widehat{\beta }$ is the vector of estimated parameters, $X$ is the matrix of input variables, and $Y$ is the vector of observed output values.

2.2. Regression Trees Ensembles: Boosted Trees, Bagging and Random Forest

Regression trees are a subset of decision trees tailored to predict continuous outcomes. These models partition the input space into regions using binary recursive splitting, where each leaf represents a specific numerical prediction (Figure 1).

The objective at each node is to select the best split by minimizing the residual sum of squares between the predicted and actual values in the resulting regions [12,13,14].

For a given dataset containing N observations with p input variables, denoted as $\left(x_i,y_i\right)$, where $x_i$ represents an input vector and $y_i$ the corresponding output, the tree construction aims to minimize the following expression (3):

$$\underset{j,s}{\mathit{min}}\left[\underset{c_1}{\mathit{min}}\sum _{x_i\in R_1\left(j,s\right)}{\left(y_i-c_1\right)}^2+\underset{c_2}{\mathit{min}}\sum _{x_i\in R_2\left(j,s\right)}{\left(y_i-c_2\right)}^2\right]$$

(3)

where $R_1\left(j,s\right)$ and $R_2\left(j,s\right)$ are the two regions after splitting on the j-th variable at point s, $c_1$ and $c_2$ are the mean predictions within these regions.

Once the input space is divided into M regions, the prediction for any new observation$x$ is the mean of the outcomes $y_i$ for the region $R_m$ to which $x$ belongs, given by (4):

$$f\left(x\right)=\sum _{m=1}^M{c}_{m}I\left(x\in R_m\right)$$

(4)

where $c_m$ represents the average output in region $R_m$ and $I\left(x\in R_m\right)$is an indicator function that is 1 if $x\in R_m$ and 0 otherwise.

Regression trees, while interpretable and flexible, tend to suffer from high variance, meaning that even small changes in the training data can result in significant differences in the structure and predictions of the tree. This instability makes regression trees prone to overfitting, especially when the model becomes too complex by learning specific noise in the training data.

2.2.1. Boosting Methodology

Boosting is an ensemble learning technique in which models are trained sequentially. Each successive model aims to improve the performance of the ensemble by correcting errors from previous iterations (Figure 2). Specifically, Gradient Boosting focuses on minimizing a predefined loss function by iteratively adding new models to an ensemble [12,13,14]. Each newly added regression tree optimizes its predictions based on the residuals, the differences between the actual and predicted values of the preceding models.

Mathematically, for a given set of training data $D=\{\left(x_i,y_i\right){\}}_{i=1}^l$, where $x_i\in R^n$ represents the input features and $y_i\in Ry$ are the target variable, the objective of Gradient Boosting is to minimize a differentiable loss function $L\left(y,\widehat{y}\right)$ where $\widehat{y}$ is the predicted output. The general form of the model can be expressed as (5):

$$F\left(x\right)=\sum _{i=1}^B{\alpha }_m{h}_i\left(x\right)$$

(5)

where $F\left(x\right)$ represents the final model prediction, ${\alpha }_m$ is the weight of each weak learner, and $h_m\left(x\right)$ is the prediction from the m-th weak learner, typically a regression tree.

At each iteration mmm, the new tree $h_m\left(x\right)$ is trained to minimize the residuals from the previous iteration. The residuals$r_i\left(m\right)$ at the m-th step are calculated as (6):

$$r_i^{\left(m\right)}=-\left[{\displaystyle \frac{\partial L\left(y_i,{\widehat{y}}_i^{\left(m\right)}\right)}{\partial {\widehat{y}}_i^{\left(m\right)}}}\right]$$

(6)

In the case of quadratic error (i.e., squared loss), the residuals simplify to (7):

$$r_i^{\left(m\right)}=y_i-F_{m-1}\left(x_i\right)$$

(7)

Here, $F_{m-1}\left(x_i\right)$ is the predicted value from the ensemble up to iteration m-1. The new regression tree is then fitted to the residuals in order to minimize the overall loss.

2.2.2. Bagging and Random Forest (RF) Methodology

To address limitation of individual regression trees and improve the predictive performance, ensemble methods such as bagging and RF are commonly used (Figure 3). These methods combine multiple decision trees to reduce variance and create more robust models. Bagging (Bootstrap Aggregating) works by generating multiple bootstrap samples from the original dataset. Each sample is created by randomly selecting data points with replacement, meaning some data points may be repeated while others are left out. For each bootstrap sample, a separate regression tree is trained independently. Instead of relying on a single tree, bagging reduces variance by averaging the predictions from all trees in the ensemble, thus providing a more stable and accurate prediction [15,16,17,18,19].

Mathematically, if $f_b^{\mathrm{*}}\left(x\right)$ represents the prediction from the b-th bootstrap sample, the final prediction from bagging is given by the average of all individual predictions (8):

$$f_{\mathit{bagged}}\left(x\right)={\displaystyle \frac{1}{B}}\sum _{b=1}^B{f}_b^{*}\left(x\right)$$

(8)

where B is the total number of bootstrap samples, and $f_b^{*}\left(x\right)$ is the prediction from the b-th regression tree. By averaging the predictions, bagging significantly reduces the model's variance, making it less sensitive to the idiosyncrasies of individual training datasets. This leads to improved generalization on unseen data, while still maintaining the flexibility and interpretability of individual decision trees.

RF build upon bagging by introducing an additional layer of randomness. Instead of considering all variables for each split, RF select a random subset of the variables at each node, which decorrelates the trees and further reduces variance. This process ensures that no single variable dominates the tree-building process.

2.3. Support Vector Regression (SVR)

Support Vector Regression (SVR) is an adaptation of the Support Vector Machines (SVM) algorithm, designed for continuous output prediction tasks. SVR transforms the input space into a higher-dimensional feature space using kernel functions, allowing for the linear separation of data points that may not be linearly separable in the original input space [20,21,22].

Given a dataset $(D=\{\left(x_i,y_i\right){\}}_{i=1}^l),$ where $x_i\in R^n$ and $y_i\in R$ are the corresponding continuous output values, the goal of SVR is to find a function $f\left(x\right)$ that has at most an $\mathsf{ϵ}$ deviation from the actual output $y_i$ for all training data while maintaining a balance between model complexity and prediction error. The regression function $f\left(x,w\right)$ in SVR is given by (9):

$$f\left(x,w\right)=\sum _{i=1}^N{w}_i\varphi \left(x_i\right)+b$$

(9)

where $\varphi \left(x_i\right)$ is a non-linear transformation of the input $x_i$ into a higher-dimensional space, $w$ is the weight vector, and b is the bias term. This model aims to minimize the following objective (10):

$$R\left(w\right)={\displaystyle \frac{1}{2}}|w{|}^2+C\sum _{i=1}^l{\mathcal{L}}_{ϵ}\left(y_i-f\left(x_i,w\right)\right)$$

(10)

where ${\displaystyle \frac{1}{2}}|w{|}^2$ is the regularization term that controls the model complexity, $C$is a regularization parameter that determines the trade-off between the model’s complexity and its tolerance to deviations, and $\left({\mathcal{L}}_{\mathsf{ϵ}}\right)$ is the loss function with an $\mathsf{ϵ}$-insensitive zone, defined as (11):

$${\mathcal{L}}_{ϵ}\left(y_i-f\left(x_i,w\right)\right)=\left\{\begin{array}{c}\begin{array}{c}0 if \left|y_i-f\left(x_i,w\right)\right|\le \epsilon \\ \end{array}\\ \left|y_i-f\left(x_i,w\right)\right|-\epsilon otherwise.\end{array}\right.$$

(11)

This loss function ensures that deviations within $\mathsf{ϵ}$ are ignored, focusing only on larger errors.

Minimizing $R$ is equivalent to minimizing (12):

$$R_{w,\xi ,{\xi }^{*}}={\displaystyle \frac{1}{2}}\left[{‖\mathit{w}‖}^2+C\left(\sum _{i=1}^l\xi + \sum _{i=1}^l{\xi }_{*} \right)\right],$$

(12)

where $\xi$ and ${\xi }^{*}$are the slack variables, which are shown in Figure 4.

The optimization problem in SVR is typically solved using Lagrange multipliers, transforming the problem into its dual form. The dual form allows the incorporation of kernel functions, which can compute inner products in the high-dimensional feature space without explicitly mapping the input vectors to that space. This is crucial for handling non-linear relationships within the data.

The dual optimization problem can be expressed as (13):

$$R_{w,\xi ,{\xi }^{*}}={\displaystyle \frac{1}{2}}\left[{‖\mathit{w}‖}^2+C\left(\sum _{i=1}^l\xi + \sum _{i=1}^l{\xi }_{*} \right)\right],$$

(13)

subject to the constraints (14):

$$\sum _{i=1}^l\left({\alpha }_i-{\alpha }_i^{*}\right)=0,\hspace{1em}0\le {\alpha }_i,{\alpha }_i^{*}\le C$$

(14)

where ${\mathsf{\alpha }}_{\mathrm{i}}$ and ${\mathsf{\alpha }}_{\mathrm{i}}^{\mathrm{*}}$ are Lagrange multipliers, $\mathrm{K}\left({\mathrm{x}}_{\mathrm{i}},{\mathrm{x}}_{\mathrm{j}}\right)$ is the kernel function that computes the dot product in the feature space without explicitly performing the transformation. The kernel function $\mathrm{K}\left({\mathrm{x}}_{\mathrm{i}},{\mathrm{x}}_{\mathrm{j}}\right)$ can take various forms which facilitate the mapping of input data into higher-dimensional spaces suitable for separation.

The prediction function for new data points $\mathrm{x}$ in SVR is then given by (15):

$$f\left(x\right)=\sum _{i=1}^l\left({\alpha }_i-{\alpha }_i^{*}\right)K\left(x_i,x\right)+b$$

(15)

Commonly used kernel functions in SVR include:

• Linear kernel $\mathrm{K}\left({\mathrm{x}}_{\mathrm{i}},{\mathrm{x}}_{\mathrm{j}}\right)={\mathrm{x}}_{\mathrm{i}}^{\mathrm{T}}{\mathrm{x}}_{\mathrm{j}}$
• Sigmoid kernel: $\mathrm{K}\left({\mathrm{x}}_{\mathrm{i}},{\mathrm{x}}_{\mathrm{j}}\right)=\tanh \left(\mathsf{\gamma }〈{\mathrm{x}}_{\mathrm{i}},{\mathrm{x}}_{\mathrm{j}}〉+\mathrm{r}\right),\mathsf{\gamma }>0$
• Radial Basis Function (RBF) kernel $(\mathrm{K}\left({\mathrm{x}}_{\mathrm{i}},{\mathrm{x}}_{\mathrm{j}}\right)=\exp \left(-\mathsf{\gamma }|{\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{j}}{|}^2\right))$.

The performance of the SVR model depends heavily on selecting appropriate hyper parameters, such as: regularization parameter C, ϵ, kernel-specific parameters, where for instance, γ in the RBF kernel, which controls the width of the Gaussian function etc. Tuning these hyper parameters is crucial for ensuring the model generalizes well to unseen data. This is typically achieved using grid search techniques to identify the optimal values for C, ϵ, and any kernel-specific parameters.

2.4. Gaussian Process Regression (GPR)

Gaussian Process Regression (GPR) is a Bayesian regression technique that offers a probabilistic framework for modeling the relationship between input features and continuous output variables. One of the significant advantages of GPR is its ability to quantify uncertainty in predictions [23].

In GPR, it is assumed that the underlying function mapping inputs to outputs is sampled from a Gaussian Process. A Gaussian Process is characterized as a collection of random variables where any finite number of these variables have a joint Gaussian distribution.

Consider the nonlinear regression problem (16):

$$y=f\left(x\right)+\epsilon , \epsilon ~{\rm N}(0, {\sigma }^2).$$

(16)

Here, $f\left(\cdot \right)$ is an unknown function mapping the input $x$ to the output $y$, and $ϵ$ represents normally distributed noise with mean zero and variance ${\sigma }^2$. GPR assumes that the unknown function $f\left(\cdot \right)$ follows a Gaussian Process, which is specified by a mean function $\mu \left(\cdot \right)$ and a covariance function (or kernel) $k\left(\cdot ,\cdot \right).$

Given a dataset with $n$ observations, $\mathrm{y}= \left\{y_1, . . . ,y_n\right\}$, these observations are modeled as a sample from a multivariate Gaussian distribution (17):

$${\left(y_1, . . . ,y_n\right)}^T~N(\mathit{\mu }, {\rm K}),$$

(17)

where:

${\mathit{\mu }=\left({\mu (x}_1\right), . . . ,{\mu (x}_n\left)\right)}^T$ is the vector of mean values.

K is the covariance matrix, where the element $K_{ij}=k\left(x_i, x_j\right)+{\sigma }^2{\delta }_{ij}$, with ${\delta }_{ij}$ being the Kronecker delta function.

For a new test input $x^{*}$, GPR predicts its corresponding output $y^{*}$. The joint distribution of both the training data outputs $y$ and the new output $y^{*}$ is a multivariate Gaussian distribution (18):

$$\left(y_1, . . . ,y_n, y^{*}\right)~N({\mathit{\mu }}^{\mathit{*}},\sum )$$

(18)

where:

• ${{\mathit{\mu }}^{\mathbf{*}}=\left({\mu (x}_1\right),...,{\mu (x}_n),\mu (x^{\mathrm{*}}\left)\right)}^T$ is the mean vector.
• $\sum$ is the covariance matrix, which can be divided into blocks (19):

  $$\sum =\left[\begin{array}{cc}K& K^{*}\\ { K}^{*T}& K^{**}\end{array}\right]$$

  (19)

In this matrix:

• $K^{\mathrm{*}}={(K\left({x^{\mathrm{*}},x}_1\right),K\left({x^{\mathrm{*}},x}_2\right),\dots ,K({x^{\mathrm{*}},x}_n\left)\right)}^T$ is the covariance between the test point and training points.
• $K^{\mathrm{*}\mathrm{*}}=K\left(x^{\mathrm{*}},x^{\mathrm{*}}\right)$ is the variance at the test point.

The conditional distribution of $y^{*}$ given the training data $y{=\left(y_1, . . . ,y_n\right)}^T$is also Gaussian, with mean and variance values, respectively (20), (21):

$${\widehat{\mathit{y}}}^{\mathit{*}}=\mu \left(x^{*}\right)+{ K}^{*T}{ K}^{-1}\left(y-\mu \right),$$

(20)

$${{\widehat{\sigma }}^{*}}^2={ K}^{**}+{\sigma }^2-{ K}^{*T}{ K}^{-1}{K}^{*}.$$

(21)

The covariance function (kernel) plays a key role in defining the properties of the functions that the GP can model. A commonly used kernel is the Squared Exponential (SE) kernel with Automatic Relevance Determination (ARD) (22):

$$k\left(x_p, x_q\right)=v^{2}exp\left[-{\displaystyle \frac{1}{2}}\sum _{i=1}^n{\left({\displaystyle \frac{x_p^i-x_q^i}{r_i}}\right)}^2\right].$$

(22)

In this equation $r_i$ controls the length scale for each input dimension, indicating how much influence each dimension has on the prediction. A large $r_i$ means the corresponding input dimension is less important.

The parameters $\left\{v,r_1, . . . ,r_n\right\}$ and the noise variance ${\sigma }^2$are called hyper parameters. These can be optimized by maximizing the log marginal likelihood of the observed data (23):

$$L\left(v,r_1, . . . ,r_n,{\sigma }^2\right)=-{\displaystyle \frac{1}{2}}\mathit{log}det K-{\displaystyle \frac{1}{2}}{ y}^T{ K}^{-1}y-{\displaystyle \frac{n}{2}}\mathit{log}2\pi .$$

(23)

2.5. Multi Gene Genetic Programming

Genetic Programming (GP) is an evolutionary algorithm-based methodology used for modeling and optimization. It generates models by simulating the process of natural selection, where a population of potential solutions evolves over time, with the best solutions "surviving" and "reproducing." This methodology is particularly suited to problems where the functional form of the model is not known beforehand.

In regression tasks, GP evolves models that predict continuous outputs based on input variables. The objective is to find a symbolic representation (a mathematical formula) that best approximates the relationship between input and output data. The flexibility of GP allows it to discover complex, non-linear relationships that may be difficult to capture using traditional regression techniques [24,25].

The process starts by generating an initial population of random solutions (Figure 5). Each individual in this population represents a possible model for the problem and it is structured using one or more decision trees. In the initial iteration, the tree model is generated by randomly selecting mathematical functions, constants, and model variables. Within this framework, each tree can be thought of as representing a single "gene" in the model. The trees culminate in terminal nodes, which are assigned either as input variables from the model or as fixed constants. All other intermediate nodes in the tree serve as functional nodes.

Each model in the population is evaluated based on a fitness function, which quantifies how well the model performs in predicting the desired output. Commonly used metrics include the Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and other error-based criteria.

Selection process of models is based on their fitness and models with better performance have a higher probability of being selected, but some randomness is introduced to maintain diversity within the population.

Selected models are paired to create new "offspring" models by exchanging parts of their structures (Figure 6). In constructing new models, the best-performing models from the previous generation are utilized to generate a new population through processes such as crossover, mutation, and direct copying. In practice, the selection of models for reproduction is probabilistic, based on their fitness score or the complexity of the model structure. Model complexity is often determined by the number of nodes and subtrees within a tree, reflecting the expressive capacity of the model.

During the crossover process, genetic material from two parent models can be exchanged entirely or partially (Figure 6). For example, suppose the model $J_1$ contains the genes $\left[G_{\left(\mathrm{1,1}\right)}, G_{\left(\mathrm{2,1}\right)}, G_{\left(\mathrm{3,1}\right)}, G_{\left(\mathrm{4,1}\right)}, G_{\left(\mathrm{5,1}\right)}\right]$ and another model $J_2$ contains the genes $\left[G_{\left(\mathrm{1,2}\right)}, G_{\left(\mathrm{2,2}\right)}, G_{\left(\mathrm{3,2}\right)}, G_{\left(\mathrm{4,2}\right)}, G_{\left(\mathrm{5,2}\right)}\right]$. If we select a random set of genes for crossover, the models may look like (24), (25):

$$J1: \left[G_{\left(\mathrm{1,1}\right)}, <{\mathit{G}}_{\left(\mathbf{2,1}\right)},{\mathit{G}}_{\left(\mathbf{3,1}\right)}> G_{\left(\mathrm{4,1}\right)}, G_{\left(\mathrm{5,1}\right)}\right]$$

(24)

$$J2: \left[G_{\left(\mathrm{1,2}\right)}, G_{\left(\mathrm{2,2}\right)},<{\mathit{G}}_{\left(\mathbf{3,2}\right)}>G_{\left(\mathrm{4,2}\right)}, G_{\left(\mathrm{5,2}\right)}\right]$$

(25)

The bolded genes represent those selected for crossover, and their exchange results in the following offspring (26), (27):

$$O_1: \left[G_{\left(\mathrm{1,1}\right)}, G_{\left(\mathrm{3,2}\right)}, G_{\left(\mathrm{4,1}\right)}, G_{\left(\mathrm{5,1}\right)}\right]$$

(26)

$$O_2: \left[G_{\left(\mathrm{1,2}\right)}, G_{\left(\mathrm{2,2}\right)}, G_{\left(\mathrm{2,1}\right)}, G_{\left(\mathrm{3,1}\right)} , G_{\left\{\mathrm{4,2}\right\}}, G_{\left\{\mathrm{5,2}\right\}}\right]$$

(27)

This high-level crossover process enables the creation of new models that inherit features from both parent models.

Additionally, crossover can occur at the gene level, known as low-level crossover, where only a portion of a gene is exchanged.

In such cases, only part of the tree is altered (as seen in Figures 6a and 6b). Mutation can also occur at this level, where a node within a gene is randomly selected, and a newly created subtree is inserted at that location (Figures 6c and 6d). These evolutionary procedures are repeated over multiple generations to refine the model.

The resulting model is pseudo-linear, as it represents a linear combination of non-linear components modeled as trees. Mathematically, this multigene regression model is represented as [24,25]:

$${\widehat{y}=b_0+b_1{t}_1+b_{2}t}_2+\dots +b_G{t}_G$$

(28)

where $b_0$ is the bias term, $b_i$ are the scaling parameters, and $t_i$are the outputs from the i-th tree (gene).

The overall gene response matrix is denoted by:

$$G=[1,t_1,t_2,\dots t_G]$$

(29)

with dimensionality $N\times (G+1))$, and the coefficient vector is:

$${b=\left[b_0{b}_1{b_2\dots b}_G\right]}^T$$

(30)

The vector b is estimated using the least-squares method, and its solution is computed as:

$$b={\left(G^{T}G\right)}^{-1}{G}^{T}y$$

(31)

Thus, the final multigene regression model can be written as:

$$\widehat{y}=Gb$$

(32)

The process of selection, crossover, and mutation is repeated over many generations. Over time, the population evolves, and the models typically improve in performance. The algorithm terminates when a stopping criterion is met. This criterion could be a pre-defined number of generations, reaching a certain fitness threshold, or convergence in the population.

2.6. Artificial Neural Networks

A Multi-Layer Perceptron (MLP) is a type of feedforward neural network designed for supervised learning tasks (Figure 7). Its architecture consists of an input layer, one or more hidden layers, and an output layer. Each layer is composed of neurons that receive signals, process them through weighted synaptic connections, and apply a non-linear activation function to generate an output [14,26].

The input layer receives the data, which is then propagated through the network. Each neuron in the hidden layers processes the inputs by computing a weighted sum of its inputs, adding a bias term, and applying a non-linear activation function such as the hyperbolic tangent function. The output layer neurons typically use a linear activation function when the model is applied to regression tasks.

Mathematically, the activation of a neuron $y_k$ a hidden layer can be expressed as (33):

$$y_k= \varphi \left.\left( \sum _{j=1}^m{w}_{kj}{x}_j+b_k\right.\right)$$

(33)

where $x_1, x_2, . . . , x_m$ are the input variables,$w_{k1}, w_{k2}, . . . , w_{km}$ represent the weights for neuron $k$, $b_k$ is the bias term, and ϕ is the activation function.

The number of neurons in the hidden layers is a critical parameter in MLP design. The selection of this number can influence the model’s performance and its ability to generalize to unseen data. The number of neurons is typically determined through experimental tuning, following various heuristic guidelines that relate the number of neurons to the dimensions of the input and output data.

The maximum number of neurons was selected as the upper limit, and various network architectures were explored, beginning with a single neuron in the hidden layer. Subsequent architectures incrementally increased the number of neurons up to the maximum value determined by the following equation (34) [12]:

$$N_H\le min\left(2N_I+1, {\displaystyle \frac{N_S}{N_I+1}} \right)$$

(34)

The choice of activation function plays a pivotal role in the learning capacity of the MLP. In this model, a hyperbolic tangent (tanh) activation function is employed in the hidden layer neurons. For the output layer, a linear activation function is utilized, which is suitable for regression problems.

The MLP model can be trained using LM algorithm, where weights are updated iteratively through gradient descent. The loss function, the mean squared error (MSE), is minimized during training to improve the accuracy of the model.

3. Dataset

The database consists of 855 test results collected from 38 research studies (Table 1), focusing on the key parameters that influence bond strength [2]. The single-lap shear test (Figure 8), recognized as one of the most effective methods for evaluating the bond strength at the FRP-concrete interface, forms the foundation of this study.

The collected dataset includes critical material properties and geometric parameters. Material properties encompass the compressive strength of a concrete cylinder (${\mathrm{f}}_{\mathrm{c}}^{\mathrm{\text{'}}}$), and the elastic modulus of FRP sheets. Geometric parameters include the thickness (${\mathrm{t}}_{\mathrm{f}}$) and width (${\mathrm{b}}_{\mathrm{f}}$) of the FRP sheet, the bond length (${\mathrm{L}}_{\mathrm{f}}$) of the FRP sheet, and the width of the concrete substrate (${\mathrm{b}}_{\mathrm{c}}$). The compressive strength of concrete cubes samples (${\mathrm{f}}_{\mathrm{c}}$) was converted to cylinder strength (${\mathrm{f}}_{\mathrm{c}}^{\mathrm{\text{'}}}$).

The filtered dataset was divided into a training set (70% of the data) for model development and a test set (30% of the data) for assessing model performance. Depending on the model, MSE (Mean Squared Error) or RMSE (Root Mean Squared Error) was used as the primary metric during training. Model accuracy on the test set was evaluated using four metrics: two absolute measures (RMSE and MAE) and two relative measures (MAPE and the correlation coefficient R). This process is illustrated in Figure 9.

The Mean Squared Error (MSE) is calculated using (35):

$$MSE={\displaystyle \frac{1}{N}}\sum _{k=1}^N{\left(d_k-o_k\right)}^2$$

(35)

where:

$d_k$ - actual value (target),

$o_k$ - predicted value,

N – total number of training samples.

The Root Mean Squared Error (RMSE) is a measure of overall model accuracy, expressed in the same units as the predicted variable (36):

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}}\sum _{k=1}^N{\left(d_k-o_{k }\right)}^2}$$

(36)

RMSE is sensitive to larger errors due to the squaring process, making it effective for identifying significant prediction discrepancies.

The Mean Absolute Error (MAE) measures the average absolute difference between predicted and actual values (37):

$$MAE={\displaystyle \frac{1}{N}}\sum _{k=1}^N\left|d_k-o_{k }\right|$$

(37)

Unlike RMSE, MAE gives equal importance to all errors, regardless of their size.

The Mean Absolute Percentage Error (MAPE), sometimes referred to as the Mean Absolute Percentage Deviation (MAPD), is a metric for assessing prediction accuracy and is calculated as follows (38):

$$MAPE={\displaystyle \frac{100}{N}}\sum _{k=1}^N\left|{\displaystyle \frac{d_k-o_{k }}{d_k}}\right|$$

(38)

In this formula (38), the absolute difference between each observed value ($d_k$) and predicted value ($o_{k }$) is divided by the observed value. This ratio is summed for all samples, averaged by dividing by the total number of samples, and finally multiplied by 100. MAPE is not suitable if any observed values are zero.

The Pearson correlation coefficient (R) is used as a relative measure to evaluate the relationship between observed and predicted values (39):

$$R={\displaystyle \frac{\sum _{k=1}^N(d_k-\overline{d}\left)\right(o_k-\overline{o})}{\sqrt{\left[\sum _{k=1}^N{\left(d_k-\overline{d}\right)}^2{\left(o_k-\overline{o}\right)}^2\right]}}}$$

(39)

where:

$\overline{o}$ - average of predicted values,

$\overline{d}$ - average of actual values.

The division of the complete dataset for machine learning models should be conducted in a manner that ensures similar statistical characteristics across subsets. This approach aims to create training, validation, and test sets with comparable distributions and statistical indicators, such as mean, standard deviation, and data range, thereby preserving data integrity and enhancing the model's ability to generalize effectively (Table 2, Table 3).

The following flowchart presents a step-by-step procedure for building a machine learning model to predict the output variable $P_u$ (kN). Each model is evaluated to determine the one that best meets the desired accuracy, making it a candidate for potential deployment.

4. Results and Discussion

A linear regression framework, enhanced not only with main effects but also with various interaction and polynomial terms, was employed to capture the complexities in bond strength prediction.

Interaction terms, such as $E_f\cdot t_f{, b}_c\cdot b_f$, $t_f{b}_f$ etc. allow us to examine how pairs of variables jointly influence bond strength, capturing multi-variable dependencies essential when one variable’s effect is modified by another.

Polynomial terms (squares of the variables), like $b_c^2$, $f_c^{\text{'}2}$, $E_f^2$ etc. addressed nonlinear relationships between predictors and bond strength, enabling the model to fit complex trends while retaining the interpretability of linear regression.

By broadening the model to include these interaction and polynomial terms, we achieved improved predictive accuracy, as validated by RMSE, MAE, MAPE, and R metrics.

The Table 4. presents the estimated coefficients, standard errors, t-statistics, and p-values for each parameter in the developed linear regression model. The most significant predictors of bond strength are the interaction between FRP thickness and sheet width ($t_f{b}_f$), the interaction between FRP elastic modulus and thickness ($E_f{t}_f$), FRP elastic modulus ($E_f$), and concrete compressive strength ($f_c^{\mathrm{\text{'}}}$), along with their interactions.

Polynomial terms indicate that the relationship between certain variables and bond strength is not strictly linear, with diminishing effects observed at higher values.

Performance metrics including Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), Mean Absolute Percentage Error (MAPE), and the correlation coefficient (R) were calculated to assess the model's (defined by equation (40)) accuracy.

The obtained regression equation is defined by the following expression (40):

$$\begin{array}{c}{P}_u= -15.92 - 0.0393 b_c+ 0.3682 f_c^{\text{'}}+ 0.0695 E_f+ 0.0654 b_f+ 0.0574 L_f+ 0.0443E_f t_f\hfill \\ + 0.0014b_c{b}_f+ 0.1313 t_f{b}_f- 0.0005f_c^{\text{'}}{L}_f- 0.0001b_c^2- 0.0025f_c^{\text{'}2}- 0.0001E_f^2\hfill \\ - 0.0009b_f^2- 0.00007L_f^2\hfill \end{array}$$

(40)

In this research, a Gradient Boosted Trees (GBT) ensemble was employed, with a grid search approach used to systematically explore and optimize the model’s key parameters. The grid search was applied to various combinations of learning rates and tree complexities to identify the best-performing configuration. Following values were investigated:

• Number of Generated Trees (NumLearningCycles = 100): The ensemble was limited to 100 trees to balance model complexity and prevent overfitting. This parameter was kept constant during the grid search to evaluate the impact of the learning rate and tree depth more precisely.
• Learning Rate (λ): Grid search was used to explore learning rates ranging from 0.001 to 1.0. The search revealed that a learning rate of 0.1 provided the best trade-off between fast convergence and error minimization.
• Number of Splits (MaxNumSplits): Tree depth, represented by the maximum number of splits, was also varied in the grid search. The search evaluated depths ranging from 1 split (shallow trees) to 512 splits (deep trees). Tree complexity was controlled by adjusting the maximum number of splits (MaxNumSplits), calculated in relation to dataset size. The maximum depth of the trees was determined using the formula ${log}_2(\mathrm{n}-1)$, where n is the number of data points and rounded to a whole number. The term n - 1 is used because the number of possible splits equals n - 1, which considers the total number of internal nodes required to split between adjacent points. Taking the logarithm of n−1 for base 2 helps determine the approximate number of splits required for full separation of the dataset. This value is then rounded to the nearest whole number to represent the maximum depth of the decision tree in terms of the number of splits. This depth calculation ensures that the trees are not too deep relative to the dataset size, which helps prevent overfitting.

The grid search identified the combination of 64 splits and a learning rate of 0.25 as the best, achieving the lowest RMSE of 3.3427 and the highest R value of 0.9536. The best MAE of 2.2542 was achieved with 128 splits and a learning rate of 0.1. The lowest MAPE was observed with 64 splits and a learning rate of 0.1, yielding a value of 0.1503 (Figure 10).

The optimal number of splits can be adopted to be 64, where the model achieved its best performance across several metrics. Beyond 64 splits, further increases in tree depth yielded diminishing returns, as seen by the leveling off of improvements in RMSE and R-squared value. Learning curve for optimal model is represented with purple color in Figure 11.

In this study, permutation error was employed to determine the relative importance of each input variable in the GBT model developed to predict the output variable $P_u$. The permutation importance method quantifies variable significance by assessing the increase in prediction error when the values of each variable are randomly shuffled, thereby disrupting its relationship with the target variable. This technique is particularly advantageous for understanding the contribution of each predictor in complex, nonlinear models like gradient boosting. The results, illustrated in Figure 12, reveal that variables $b_f$ (width of the FRP) and $L_f$(bond length) have the highest importance, significantly influencing the model's accuracy. In contrast, variables such as $b_c$ (concrete width) and $f_c^{\text{'}}$ (concrete compressive strength) have a lesser impact on the predictive power.

Bagging (Bootstrap Aggregation) is the primary method used to implement the ensemble model in this research. Depending on whether we use only a subset of the input variables or the entire set of variables, a RF or TreeBagger model will be created, respectively. This method involves creating multiple decision trees (500 trees in your case) on different bootstrap samples of the data and averaging their predictions to improve model robustness and reduce variance.

The research specified 500 learning cycles (NumLearningCycles = 500), meaning 500 base decision trees are trained, and their predictions are aggregated to form the final output. This large number of trees stabilizes the model and helps capture complex patterns in the data.

Two hyper parameters are tuned during the training process: Min Leaf Size and Number of Variables to Sample at each split.

• Min Leaf Size: This parameter controls the minimum number of observations required to form a leaf in a decision tree. Smaller leaf sizes result in deeper trees, allowing the model to capture more detailed patterns in the data but also increasing the risk of overfitting. In this implementation, Min Leaf Size is varied from 1 to 10.
• Number of Variables to Sample: At each split in a decision tree, a subset of predictor variables is randomly selected for consideration. The number of variables sampled at each split (NumVariablesToSample) is varied from 1 to 6.

Across all metrics (RMSE, MAE, MAPE, and R), smaller Min Leaf Size values (especially 1 and 2) consistently lead to better results (Figure 13). This suggests that deeper trees, which allow for more detailed splits, are necessary to capture the complexity in the data. The best-performing models used 3 or 4 variables for splitting at each node. This aligns with the idea that a moderate number of variables strikes a balance between simplicity and model complexity.

The best configurations identified are:

• Best RMSE: 3.7860, with Min Leaf Size = 1 and Number of Variables = 3.
• Best MAE: 2.5821, with Min Leaf Size = 1 and Number of Variables = 4.
• Best MAPE: 0.1778, with Min Leaf Size = 1 and Number of Variables = 4.
• Best R-squared: 0.9415, with Min Leaf Size = 1 and Number of Variables = 3.

There is consistency in the best-performing configurations across the different evaluation metrics (RMSE, MAE, MAPE, and R). The model with Min Leaf Size = 1 and Number of Variables = 3 or 4 appears to be a robust choice for optimal performance across various metrics.

The SVR used a radial basis function (RBF), linear and sigmoid kernel and applied grid search cross-validation to optimize the hyper parameters of the SVR model. All data were normalized. Normalization scales the data to a range between 0 and 1 by subtracting the minimum value and dividing by the range (max-min). This step is critical in SVR because the optimization process is sensitive to the scale of the input data. After normalization, the training and testing data were used in the SVR model.

Grid search is used to find the optimal values for three SVR hyper parameters:

• C (cost/regularization): This controls the trade-off between allowing slack variables (errors) and forcing the decision boundary to be as tight as possible. A higher C makes the model focus more on correctly classifying all training points but risks overfitting.
• Gamma (γ): This defines the influence of individual training examples. Smaller values of gamma imply that each training point has a far-reaching influence, while higher values imply more localized influence.
• Epsilon (ε): This defines a margin of tolerance where no penalty is given to errors within a certain range. Epsilon controls the sensitivity of the model to prediction errors.

First, the analysis of the application of the RBF kernel was performed. The parameters ${log}_2\left(C\right)$, ${log}_2\left(\mathsf{\gamma }\right)$, and ${log}_2\left(\mathsf{\epsilon }\right)$ are varied from -5 to 5 in coarse increments of 1. For each combination of C, $\mathsf{\gamma }$, and $\mathsf{\epsilon }$, an SVR model is trained, and the MSE is calculated on the test set. The best combination of hyper parameters is identified by minimizing the MSE.

Finer grid search is applied around the best values of C, gamma, and epsilon found in the coarse grid search. The step size is reduced to half the previous value, and a smaller search range is applied around the best initial values. The model is iteratively again retrained for each combination, and MSE is recalculated to find the refined best parameters. The accuracy criterion values for the analyzed models are given in Table 5.

An analogous procedure was applied to the other two kernel functions. The optimal models were determined with following kernel parameters:

• C = 1.4513; ε = 0.0043; γ = 20.7363 for the RBF kernel;
• C = 0.3208 and ε = 0.0432 for the linear kernel;
• C = 23.6326; ε = 0.0521; γ = 0.0118 for sigmoid kernel.

In this research, various covariance functions were employed to design the GPR model, each with different strengths in modeling the relationships between the input variables and the target. The covariance functions explored included the exponential, squared exponential, Matérn (with 3/2 and 5/2 degrees of smoothness), and rational quadratic kernels. These kernels control how the model captures correlations in the data, effectively shaping how predictions are influenced by nearby points. Additionally, Automatic Relevance Determination (ARD) versions of these kernels were tested, which assign individual length scales to each input variable. This feature allows the model to automatically determine the relative importance of each predictor, improving performance in datasets with features of varying significance.

To ensure comparability between features and prevent any single feature from disproportionately affecting the model, the input data was standardized. The standardization was achieved using Z-score normalization, where each feature is rescaled to have a mean of zero and a standard deviation of one.

This transformation is essential in GPR models to ensure that all variables contribute appropriately to the prediction process, especially when using covariance functions sensitive to the scale of the input data.

The GPR model employed a constant basis function, meaning the model assumes a fixed baseline for the predictions, allowing the kernel functions to capture the complexity in the relationships. The model was trained and evaluated using a holdout validation approach.

The parameters of the covariance functions used in the analysis were optimized through gradient-based methods, applied to the log marginal likelihood expression (Table 6, Table 7)). The length scale values for each input variable offer insights into the importance of these variables concerning the model’s predictive performance (Table 7). Variables with smaller length scales tend to have a greater influence on the model, as the relationship between the input and output changes more rapidly with those features. In contrast, larger length scales indicate less significant features.

All models with ARD perform better than their non-ARD counterparts (Table 8, Table 9). The ARD mechanism allows the models to focus on the most important features, leading to better predictive performance and lower error metrics.

The provided Figure 14 uses the inverse of these length scales to visually highlight variable importance, the smaller the length scale, the larger the corresponding bar in the graph, indicating greater influence on the output. The results indicate that among all the input variables, the bond length of FRP ($L_f$) plays the most crucial role in predicting bond strength. This could be due to the strong relationship between the bond area and the bond strength. The thickness of FRP ($t_f$) and concrete compressive strength ($f_c^{\text{'}}$) are also important but to a lesser extent, aligning with structural engineering principles where material properties impact bonding but not as strongly as the bond length. The relatively low importance of $E$ and $b_c$ suggests that variations in these factors do not significantly alter the model’s predictions for bond strength, potentially because their influence is overshadowed by the more dominant variables.

In this study, it is implemented a Multigene Genetic Programming (MGGP) approach that involved a comprehensive analysis of different values for the number of genes and varying tree depths, both of which are critical parameters that determine the complexity and performance of the resulting models. The method aimed to explore a wide range of configurations to identify the optimal balance between model complexity and predictive accuracy.

The MGGP method was designed to iteratively search through the solution space by repeating the procedure ten times in this research. This iterative repetition was necessary to address the inherent randomness in the initialization of model parameters. Each iteration began with a random set of parameters, resulting in a unique model at the end of each run.

The parameter settings used in this study are summarized in Table 10, which provides key details such as the function set, population size, number of generations, maximum number of genes, and maximum tree depth. These parameters were carefully chosen to enable an extensive exploration of possible model structures while maintaining control over computational complexity. Specifically, models were analyzed with the number of genes ranging from one to six and tree depths varying from one to six. These settings allowed the study to capture a wide range of model structures, from simple to highly complex.

During the training process, two main objectives were prioritized: minimizing RMSE and minimizing the expression complexity of the models. These two factors were treated as part of a multi-objective function, balancing prediction accuracy with interpretability. RMSE was used as a standard measure with weight of 0.7 of model accuracy, indicating how well the model predictions aligned with the actual data. Meanwhile, expression complexity with weight of 0.3 was an important metric for ensuring that the resulting models remained interpretable, avoiding overly complex expressions that would be difficult to understand or deploy.

To enhance the robustness of the results, the process was repeated over 1000 generations, allowing the evolutionary process to refine the models over time. The population size initially was set to 100, ensuring that a diverse range of solutions could be explored in each generation. After all number of generations is achieved process is restarted and repeated and new models are created. These multiple models, derived from ten independent runs, were ultimately merged to form one final population that was representative of the most promising solutions.

These models were then subjected to a detailed evaluation, focusing on their correlation coefficients and expression complexities. The models that exhibited superior correlation coefficients, indicating strong predictive power, and lower complexity were selected for further analysis. This selection formed what is known as the Pareto front (as depicted in Figure 11). The Pareto front represents the set of models that provide the best possible trade-offs between competing objectives namely, accuracy and complexity. By selecting models along this front, we ensured that the final set of models included those that were both effective in making accurate predictions and interpretable due to their reduced complexity.

The use of a tournament selection mechanism, with a size of 2, allowed for a balance between exploration and exploitation of the solution space. The elitism rate of 0.05% ensured that the best individuals from each generation were preserved, preventing the loss of high-quality solutions. Crossover and mutation probabilities were set to 0.84 and 0.14, respectively, to encourage diversity within the population and prevent premature convergence.

The MGGP approach used in this study involved an extensive analysis of different model configurations, focusing on balancing accuracy and complexity. By iteratively generating and evaluating models, and ultimately selecting those along the Pareto front, we were able to identify models that were both accurate and interpretable (Figure 15). This approach ensured that the final model was not only effective in terms of predictive performance but also practical for real-world applications due to its manageable complexity.

The analytical expression of the optimal Pareto front model (Model ID 7961) is defined by the following equation (41):

$${\sigma =b+\sigma }_1+{\sigma }_2+{\sigma }_3+{\sigma }_4=0.53{ + \sigma }_1+{\sigma }_2+{\sigma }_3+{\sigma }_4+{\sigma }_5+{\sigma }_6$$

(41)

Analytical expressions for individual genes are given in the Table 11.

Each of the symbolic analytical expressions or genes listed in the Table 11 is also structurally represented in the form of a tree in Figure 16.

The relationship between the accuracy of the obtained model which consists of a population of 800 individuals and two parameters, namely the Number of Genes and the Depth of the Tree, is illustrated in the Figure 17. The Pareto front depicted in the Figure 15 is defined using two criteria: the values RMSE and the coefficient of determination R². It represents a set of non-inferior (Pareto-optimal) solutions based on these criteria. The GPTIPS software used for model creation in this research highlights the optimal model according to RMSE and R² criteria with a green circle outlined in red. Since the analysis in this study was conducted using the criteria of RMSE, MAE, MAPE, and Correlation Coefficient R, each model from the aforementioned set of non-inferior solutions was further analyzed concerning all defined criteria. This was done to identify a potential model with lower complexity.

From the Table 12, it can be observed that the optimal model from the set of non-inferior solutions, with a complexity of 95, has very similar accuracy to the model labeled as Model ID 6959, which has a complexity of 42. Regarding the three defined criteria MAE, MAPE, and R, the differences are almost negligible, while in terms of RMSE, the difference appears in the first decimal place.

By applying this lower-complexity model, a mathematical expression significantly simpler for practical use is obtained, whose accuracy is only slightly less than that of the optimal model.

By employing the MGGP model using the GPTIPS software [22 ,65}, models were developed where only variables significant to the problem which is researched were retained through evolutionary processes.

Throughout this evolutionary process, variables that do not significantly contribute to the model's accuracy are generally eliminated, resulting in a final model comprising only the most relevant predictors. However, in this particular analysis, all six input variables were found to contribute to the model's accuracy and thus need to be incorporated into the model formulation.

The analytical expression of the model (Model 6959) with reduced complexity is defined by the following equation (42):

$${\sigma =b+\sigma }_1+{\sigma }_2+{\sigma }_3+{\sigma }_4=0.53{ + \sigma }_1+{\sigma }_2+{\sigma }_3+{\sigma }_4+{\sigma }_5+{\sigma }_6$$

(42)

Analytical expressions for individual genes for Model 6959 (MGGP simplified) are given in Table 13.

The neural network model architecture was established through a trial-and-error method. In this study, the input layer comprises 6 neurons, and the output layer consists of a single neuron.

The optimal number of neurons in the hidden layer was determined experimentally. Beginning with a single neuron, the number of neurons was incrementally increased, with each configuration assessed based on RMSE, MAE, R, and MAPE metrics. Models utilizing the Levenberg-Marquardt (LM) algorithm were tested, with data divided into same training and test sets for analysis.

All model variables underwent linear transformation to ensure all variables were scaled to the range [−1,1]. Here, the minimum value was mapped to -1 and the maximum to 1, with intermediate values scaled linearly. Scaling was applied to equalize the influence of each variable, as a variable’s absolute size may not reflect its actual impact. The same standard settings in MATLAB were applied across all model architectures during training (Table 14).

As guidelines for estimating the upper limit of the neuron count in the hidden layer, the previously mentioned expressions were used, with a recommendation to adopt the lower value as follows:

$$N_H\le 2{\times N}_i+1=2\times 6+1=13$$

$$N_H\le {\displaystyle \frac{N_s}{N_i+1}}={\displaystyle \frac{855}{6+1}}=122.14.$$

The RMSE, MAE, MAPE, and R values for different network architectures, tested with up to 23 neurons, are illustrated in the Figure 18.

These metrics indicate that the model with 13 neurons in the hidden layer achieved optimal performance, close to the upper boundary of the initially tested range. Consequently, the range was expanded to evaluate the impact of additional neurons. However, adding more neurons did not enhance accuracy; rather, it highlighted a risk of overfitting. Beyond this limit, the model began to fit too closely to the training data, reducing its ability to generalize effectively to unseen data.

A comparative analysis of all analyzed models according to the adopted accuracy criteria is given in the Table 15.

Table 15 presents a comparative analysis of several machine learning models used for predicting bond strength, evaluated based on four key accuracy metrics: Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), Mean Absolute Percentage Error (MAPE), and the correlation coefficient (R). Among the models evaluated, the GPR ARD Exponential model outperforms all others, with the lowest RMSE (2.8671 MPa), MAE (1.9319 MPa), and MAPE (0.1329%) while achieving the highest correlation coefficient (R = 0.9658). This indicates its superior accuracy and predictive power.

The Gradient Boosted model also shows strong performance, with an RMSE of 3.3427 MPa and an R of 0.9536, making it a competitive alternative for situations requiring complex, non-linear modeling. Similarly, RF and TreeBagger models also perform well, with RMSE values of 3.7860 MPa and 3.8302 MPa, respectively, and R values above 0.939. These models are effective in reducing variance and offer robustness through ensemble learning.

The SVR RBF model, while slightly less accurate than the tree-based models, still achieves reasonable results with an R of 0.9332 and a relatively low MAE (3.5352 MPa), suggesting its applicability in scenarios where a simpler model is preferred.

On the other hand, the MGGP and MGGP simplified models, despite their interpretability, show lower accuracy with higher RMSE values (4.4436 MPa and 4.6829 MPa, respectively). This indicates that while these models may be useful for applications requiring simple, symbolic expressions, they may not be the best choice when high accuracy is critical.

Finally, the Neural Network (NN 6-13-1) model, with an RMSE of 4.0992 MPa and R of 0.9293, offers a balanced trade-off between model complexity and accuracy. However, it still falls short compared to the top-performing models like GPR and Gradient Boosted.

In summary, the GPR ARD Exponential model is the optimal choice based on all evaluated metrics, though simpler models such as Gradient Boosted or RF may still be valuable depending on the application, especially in terms of balancing complexity and performance.

Figure 19 shows a regression plot, comparing the target values with the predicted values from the model. This plot illustrates the accuracy of the model's predictions, with the closer the predicted values are to the target values, the better the model's performance. The regression plot reflects a strong correlation between the actual and predicted results, indicating the model's ability to effectively capture the relationship between the input features and bond strength.

Figure 20 provides a direct comparison between the target and predicted values for the optimal GPR ARD Exponential model. This visual further demonstrates how well the model predicts the bond strength, with the predictions closely matching the target values across the dataset. The closeness of the predicted values to the target values confirms the high accuracy and reliability of the model in predicting bond strength.

While the model is relatively complex, a MATLAB Graphical User Interface (GUI) was developed to facilitate practical implementation by experts. The model's code includes a built-in database, making it easy to adapt the model through the expansion of the database with additional data. The MATLAB code is provided in the Appendix A1. The GUI functions by allowing the user to input the corresponding values for the six input variables, and the prediction is obtained by clicking the 'Predict Bond Strength' button (Figure 21).

5. Conclusions

In this research, a comprehensive evaluation of machine learning models was conducted to predict the interfacial bond strength of fiber-reinforced polymer (FRP) applied to concrete. Using an extensive dataset of 855 experimental results from 38 studies, the study assessed multiple advanced algorithms, including multiple linear regression with interactions, ensemble methods (Gradient Boosted Trees, TreeBagger, Random Forest), Support Vector Regression (SVR), Gaussian Process Regression (GPR), Multigene Genetic Programming (MGGP), and neural networks.

Among the models evaluated, GPR with Automatic Relevance Determination (ARD) proved to be the most accurate, achieving an RMSE of 2.8671 MPa, MAE od 1.8953 MPa and a correlation coefficient (R) of 0.9658. This model identified the FRP bond length, FRP thickness, and concrete compressive strength as the most significant predictors of bond strength, emphasizing their roles in effective FRP-concrete reinforcement.

The Gradient Boosted Trees ensemble model also showed high accuracy with an RMSE of 3.3427 MPa, MAE of 2.2603 MPa and R of 0.9536, especially when optimized through grid search for learning rate and tree depth, making it a strong candidate for complex, non-linear interactions between variables. While MGGP achieved slightly lower accuracy, it provided interpretable models that balance predictive performance with simplicity, which is advantageous for practical applications.

This analysis demonstrated that the choice of machine learning model impacts predictive accuracy, with ensemble methods and GPR offering robust, high-performing options for predicting FRP-concrete bond strength.

Additionally, by leveraging interaction terms and polynomial features for linear regression, the study enhanced the models' ability to capture the non-linear relationships inherent to FRP-concrete interactions, resulting in a direct and straightforward expression that is practically applicable.

This research underscores the potential of machine learning to enhance the reliability of FRP applications in structural engineering, contributing valuable insights to civil engineering reinforcement techniques.