1. Introduction

The introduction of recycled aggregates into 3D concrete printing advances sustainable and automated construction by making use of recycled materials and reducing the consumption of natural resources and waste during construction. However, this poses challenges in terms of the material itself as well as the process^[1-3]. At the microscale, residual mortar adhering to the recycled aggregates’ surfaces increases porosity and weakens the interfacial transition zones. At the macroscale, the layer-by-layer deposition process creates numerous inter-layer interfaces, resulting in pronounced anisotropic behavior^[4]. This makes the prediction of the splitting tensile strength (STS) of 3D printed recycled aggregate concrete (3DPRAC) challenging due to the difficulties involved in using conventional empirical or experimental methods.

Accurate prediction of the mechanical performance in 3D printed concrete requires materials that satisfy specific rheological criteria for printability^[5]. Therefore, this study focuses on mixtures with validated printability to accurately predict their hardened mechanical properties. Nevertheless, conventional empirical models and isolated experimental approaches remain inadequate for accurate prediction, even within this rheologically qualified range. This inadequacy arises from the complex nonlinear coupling between material parameters, such as the water-to-binder (w/b) ratio and recycled aggregate replacement ratio (RCA-r), and process-related variables, including layer height and loading direction. Therefore, advanced data-driven approaches capable of capturing nonlinear relationships and providing interpretable predictions are essential to comprehensively evaluate the combined effects of material constituents and printing process parameters on the mechanical properties.

Data-driven computation is growing in popularity to tackle such intricacies while bypassing the prohibitive costs of iteratively executing physical experiments^[6-11]. Techniques from the broader soft computing domain have proven to be effective for reducing experimental workload while achieving accurate results^[12,13]. Specifically, algorithms such as Convolutional Neural Networks^[14,15] and gene expression programming (GEP)^[16,17] have been employed for modeling recycled aggregate concrete (RAC). Researchers have successfully deployed diverse models for STS prediction, including Artificial Neural Networks^[16,18-23], Support Vector Machines^[24], and ensemble methods such as random forest (RF)^[25,26] and XGBoost^[21,27]. Furthermore, recent works have pointed out the great importance of integrating automatic hyperparameter optimization for robust model development, as this mechanism eliminates manual tuning bias and accurately identifies globally optimal configurations^[28,29].

Despite its immense potential, three critical limitations hinder the application of machine learning to the complex 3DPRAC system. First, current models exhibit notable insensitivity to printing process parameters. These models focus primarily on material composition mapping while neglecting the unique characteristics of 3D concrete printing. This oversight prevents the quantitative description of the significant anisotropy induced by variables such as layer height and loading direction. Second, existing studies lack generalization capability across diverse material systems, as they are typically confined to singular material domains. Consequently, these models fail to discern universal evolutionary patterns when confronted with heterogeneous datasets comprising RAC, fiber-reinforced recycled aggregate concrete (FR-RAC), and 3DPRAC, resulting in a significant reduction in prediction accuracy. Third, the “black-box” nature of algorithms is a major obstacle to their uptake. Although deep neural networks achieve acceptable accuracy, the lack of a physical or mathematical expression leads to poor interpretability, thereby severely restricting trust and widespread adoption in engineering practice.

To address these challenges, this study develops a prediction model for the STS of 3DPRAC that integrates mixture proportions and printing process parameters. GEP is employed as the modeling technique due to its capability to automatically generate explicit mathematical expressions through symbolic regression, offering both transparency and direct applicability in engineering practice. The proposed model incorporates printing layer height and loading direction as key input variables, combined with feature engineering, to enhance interpretability and reveal the underlying mechanisms governing the mechanical behavior of 3DPRAC.

The rest of this paper is structured as follows (as depicted in Figure 1). Section 2 describes the dataset collection and preprocessing, and the GEP model development. Section 3 covers the model performance evaluation and the sensitivity analysis of key parameters. Finally, Section 4 concludes the study.

Display Full Size

Figure 1. Organizational framework of the study. 3DPRAC: 3D printed recycled aggregate concrete; FR-RAC: fiber-reinforced recycled aggregate concrete; RAC: recycled aggregate concrete; STS: splitting tensile strength; GEP: gene expression programming; CGEP: chaotic gene expression programming.

2. Methodology

2.1 Data collection and description

This study aggregated a dataset of 110 data points from two sources^[30-32]. The first database consists of data derived from peer-reviewed literature on RAC and FR-RAC systems. The second database was derived from prior experimental work conducted by the authors on 3DPRAC.

This section introduces a dataset construction method for STS prediction of multiple materials. This method differs from traditional methods that rely on learning individual parameter categories. RCA-r and fiber volume (FV) are set as continuous variables to predict global variable relationships, significantly improving the predictive robustness of small sample data models and reducing the risk of overfitting.

The dataset consists of six variables, as presented in Table 1. Load direction (LD) is a binary variable used to represent the anisotropic mechanical behavior of 3D printed concrete, where a value of 0 indicates loading parallel to the printing layer, while a value of 1 indicates loading perpendicular to the printing layer, thereby denoting the parallel and perpendicular loading directions. The average value of 0.35 indicates that approximately 35% of the samples were loaded perpendicular to the printing layer. This mirrors the distribution of the actual data, where RAC and FR-RAC are labeled LD = 0 (i.e., isotropic), while 3DPRAC specimens are classified according to their tested direction. The dataset was randomly sampled and divided into a 70% training set and a 30% testing set.

2.2 Data preprocessing

2.2.1 Feature engineering

Data preprocessing and feature engineering constitute the cornerstone of predictive accuracy and interpretability in machine learning applications. A physics-informed variable processing strategy was formulated to leverage the capacity of GEP for extracting explicit mathematical expressions. Capturing of the underlying physical mechanisms governing materials and processes is achieved by converting discrete experimental descriptions into continuous input signals with distinct physical implications.

The continuous mapping of material parameters was realized by leveraging the decoupling strategy described in Section 2.1, taking RCA-r and FV as continuous numerical inputs. This yields a continuous physical response surface over a range of material systems, allowing for an analytical description of the nonlinear relationships between mixture proportions and mechanical strength.

LD is a variable reflecting the response differences of 3D printing materials to loading directions. When loads are applied parallel to the printing direction (LD = 0), the continuous material strength is primarily tested, whereas when applied perpendicular to the printing direction (LD = 1), the performance is related to interlayer properties. This assignment encodes a physical hypothesis where an LD value of 1 triggers interlayer weakening. Accordingly, the negative LD coefficient in the GEP-evolved formula captures the strength loss due to interlayer defects.

The dimension of time was introduced by making age a predictor variable, allowing for the learning of the decay rate of cement hydration using a logarithmic or exponential function. Finally, min-max normalization (Eq. (1)) was applied to scale all continuous features to the range [0, 1], reducing scale differences and allowing them to contribute equally during training, which accelerated convergence.

(1)$X_i=\frac{X_i-X_{min}}{X_{max}-X_{min}}(i=1,2,\dots ,n)$

where X_min and X_max denote the minimum and maximum values of the feature X, respectively.

2.2.2 Statistical analysis of variables

Figure 2 shows a correlation heatmap to assess inter-variable dependencies and to evaluate the degree of multicollinearity. Cell color intensity indicates the strength of correlations between variables. Layer height (LH) demonstrated a strong positive correlation with LD (r = 0.90), while the w/b ratio showed a strong negative correlation with LH (r = -0.71). Most other absolute correlation coefficients remained below 0.3, thereby avoiding potential multicollinearity issues typically associated with strong correlations (|r| > 0.5). This overall low level of correlation suggested that this dataset had a reasonable level of independence and would not introduce bias into model training.

Display Full Size

Figure 2. Correlation heatmap analysis of variables. FV: fiber volume; STS: splitting tensile strength; RCA-r: recycled aggregate replacement ratio; LH: layer height.

The linear correlations between the target variable STS and individual input variables were uniformly low, with the maximum coefficient reaching only 0.32 for Age, and others remaining below 0.2. This observation underscored the inadequacy of single physical parameters in elucidating strength variations, thereby necessitating the adoption of holistic multi-parametric analysis and non-linear modeling approaches to capture the complex synergistic effects among variables.

Variance inflation factor (VIF) and Tolerance were calculated to evaluate the severity of multicollinearity, as detailed in Table 2. All VIF indices remained significantly below the critical threshold of 10, indicating negligible multicollinearity risks, despite marginal elevations observed for w/b (3.549) and RCA-r (3.214). Tolerance values for FV and LH exceeded 0.75, supporting the robust independence among predictors.

The LD variable was deliberately excluded from the multicollinearity assessment based on two pivotal considerations. First, regarding statistical applicability, the discrete categorical nature of LD renders it incompatible with VIF metrics intended primarily for continuous linear dependencies. Second, an inherent structural coupling between LD and the material system stems from the contrast between the single-direction constraint of non-3D printed specimens and the multi-directional nature of 3D printed counterparts. This distinction generates a correlation artifact stemming from experimental design rather than genuine physical collinearity. Accordingly, LD was retained as a distinct control parameter to ensure the decoupled elucidation of anisotropic effects on STS.

Both correlation analysis and VIF diagnostics confirmed that the selected features met the statistical assumptions for multivariate regression, providing a sound basis for model development.

2.3 Gene expression programming

2.3.1 GEP working mechanism

GEP is an advanced evolutionary computation methodology based on the constructs of natural selection and genetics. Resulting candidate solutions are coded in linear chromosomes of fixed length called the genotype, which are then mapped to non-linear expression trees called the phenotype (Figure 3). GEP achieves a balance between rapid convergence and superior generalization through the unconstrained evolution of mathematical expressions that are intelligible and explicit^[33,34]. The chromosome assumes a head-tail architecture, with the head containing both function and terminal symbols and the tail containing terminal symbols only. This structure means that chromosomes that have been changed by mutation or crossover can map directly onto syntactically valid expression trees, thus allowing the free optimization of challenging mathematical regularities.

Display Full Size

Figure 3. Expression tree of chromosome.

The GEP model was created through an iterative optimization process, illustrated in Figure 4. It began with the random generation of individuals, all of which were expressed as chromosomes in Karva notation. The algorithm then followed with a core evolutionary loop whereby a fitness function assessed each phenotype, and the principle of the survival of the fittest dictated that the finest candidates selected served as the breeding pool. Then, diversity was retained in the population, and the solution space was elaborated by genetic operations of replication, crossover, and mutation until user-defined termination criteria were satisfied.

Display Full Size

Figure 4. GEP working mechanism. GEP: gene expression programming.

Notably, model performance is sensitive to the proper tuning of hyperparameters, including population size, head length, and probabilities for the genetic operators. The complex physico-mechanical behavior of 3D printed recycled concrete was captured using diverse functions, including arithmetic, algebraic, and nonlinear operators.

In this study, an STS prediction model was established with the GEP for multiple material systems, demonstrating great precision for 3D printed recycled concrete among others. It was essential to transform complex data-driven results into achievable design tools, providing explicit mathematical expressions that were easy to apply directly in engineering, thus offering strong theoretical support for mix design and performance prediction.

2.3.2 Model design and parameter configuration

This study developed two GEP models with distinct hyperparameter configurations (designated CGEP1 and CGEP2) using the GeneXproTools 5.0 platform. The objective was to systematically optimize prediction accuracy and evaluate the influence of model complexity. CGEP1 served as the baseline, featuring 4 genes and 60 chromosomes with a function set of basic arithmetic operators (such as +, -, *, and /). CGEP2 enhanced solution quality through structural adjustment by augmenting the gene count to 5, setting the population size to 50, and incorporating square root (Sqrt) and natural logarithm (Ln) functions to capture nonlinear characteristics. A head length of 12 was maintained for both models, with addition utilized as the linking function to integrate substructures derived from individual genes. Detailed hyperparameter settings are listed in Table 3.

The effective implementation of GEP hinges on the precise definition of five core components, comprising the function set, terminal set, fitness function, control parameters, and termination criteria. An initial population was randomly generated based on the selected function and terminal sets. Chromosomes encoded expression trees of varying topologies within the head and tail architecture. Sub-expression trees (Sub-ETs) from individual genes were aggregated into a final mathematical model via a linking function. Throughout the evolutionary process, root mean square error (RMSE) served as the fitness function to quantify individual performance. Subsequently, the algorithm iteratively optimized the population through genetic operators, including reproduction, crossover, and mutation, until the maximum iteration count was reached or model convergence was achieved.

2.4 Model evaluation indicators

Three standard statistical criteria were adopted for assessing the goodness of fit of the GEP models: the Coefficient of Determination (R²), mean absolute error (MAE), and RMSE. These are defined mathematically by Eqs. (2) to (4):

(2)$R^2=1-\frac{\sum _{i=1}^n{(X_i-Y_i)}^2}{\sum _{i=1}^n{(\overline{X}-Y_i)}^2}$

(3)$MAE=\frac{\sum _{i=1}^n|X_i-Y_i|}{n}$

(4)$RMSE=\sqrt{\frac{\sum _{i=1}^n{(X_i-Y_i)}^2}{n}}$

where, n denotes the total sample size, while X_i and Y_i represent the predicted and actual values for the i-th sample, respectively, with $\overline{X}$ indicating the arithmetic mean of predicted values. R² ranges from 0 to 1, with values tending towards 1 indicating better agreement between the predictions and the observed data. MAE indicates the average absolute error between predictions and observations, with lower values implying a better fit. RMSE is the square root of the mean squared error and is more sensitive to large deviations, making it useful for evaluating model performance on outliers.

Overall, R² measures linear correlation, while the MAE and RMSE quantify the predictive discrepancies. The comparison of performance outcomes across training and testing datasets effectively diagnosed potential overfitting, ensuring robust generalization ability to unknown data.

2.5 Sensitivity analysis

This study used both the relevance factor (r) and permutation feature importance (PFI) to systematically assess the relative effects of individual input variables on model predictions.

2.5.1 Relevance factor analysis

The relevance factor assesses the linear correlation strength between input variables and model outputs. It elucidates the relative statistical significance of each parameter by quantifying the direct impact of variations in individual input parameters on STS predictions. The mathematical definition is presented in Eq. (5):

(5)$r(I_j,O)=\frac{\sum _{i=1}^N(I_{j,i}-{\overline{I}}_j)(O_i-\overline{O})}{\sqrt{\sum _{i=1}^N{(I_{j,i}-{\overline{I}}_j)}^2\sum _{i=1}^N{(O_i-\overline{O})}^2}}$

where I_j,i is the j-th input value of the i-th sample and O_i is the predicted STS output for the i-th sample. ${\overline{I}}_j$ and $\overline{O}$ represent the arithmetic means of the j-th input parameter and model output, respectively, while N indicates the total sample size. Higher absolute r values signify stronger explanatory power of the input parameter with respect to model output, confirming its significance in STS prediction; conversely, lower values imply weaker influence. To ensure analytical robustness, calculations were conducted using the training dataset, adhering to methodologies defined in existing literature^[35,36].

2.5.2 Permutation feature importance

PFI evaluates feature importance by randomly permuting feature values and observing the decline in model performance. The greater the performance degradation, the more important the feature is, as this reveals its contribution to predictions. This is calculated as defined in Eq. (6):

(6)$PFI=MA{E}_p-MA{E}_o$

where MAE_o is the mean absolute error on the original dataset and MAE_p is the mean absolute error obtained by randomly permuting a specific input variable. Higher PFI values indicate stronger variable influence on the model predictions, whereas lower values indicate weaker influence. Values close to zero suggest that the corresponding input variables have little impact. PFI can effectively identify key variables that affect STS performance, providing a foundation for subsequent model refinement and feature selection.

3. Results and Discussion

This section constitutes a comprehensive evaluation of the prediction capabilities of the GEP-based unified STS framework. First, the overall prediction accuracy and error distribution are statistically analyzed and compared. Second, sensitivity analysis reveals the influence of input variables, particularly 3D printing parameters and loading direction, on STS. Finally, the explicit GEP formula is analyzed to demonstrate the model’s interpretability and effectiveness.

3.1 Overall prediction performance of GEP models

As detailed in Section 2.3.2, two GEP models featuring distinct topological architectures and hyperparameter settings, designated as CGEP1 and CGEP2, were constructed. This section aims to identify the optimal model through quantitative evaluation on the testing dataset.

3.1.1 Statistical metrics and regression analysis

R², RMSE, and MAE were used as performance metrics to evaluate predictive accuracy. Table 4 presents the performance metrics for CGEP1 and CGEP2 on the training and testing datasets.

As shown in Table 4, while both models demonstrated reasonable fitting proficiency during training, performance discrepancies on the testing set reveal fundamental differences in generalization potential. The lower R² obtained for CGEP1 suggests limited predictive capacity for unseen data, while elevated RMSE and MAE values confirm significant deviations when processing new samples. By contrast, CGEP2 exhibits superior predictive fidelity. An R² of 0.914 on the testing set indicates that this model captures over 91% of the variability in the target variable, achieving high concordance between predicted and actual values. Moreover, CGEP2 achieves significantly lower RMSE and MAE values than CGEP1. This improvement stems from the increased number of genes, extended head length, and the inclusion of non-linear operators (e.g., Sqrt), which enable CGEP2 to better capture the complex non-linear relationships between input variables and STS.

Figure 5 presents the regression analysis between predicted and experimental values for CGEP1 and CGEP2. The figure includes data from both training and testing sets, as well as the ideal diagonal (y = x) and ± 5% error bands. Figure 5a shows that CGEP1 exhibits a scattered distribution, with many data points deviating from the 1:1 line and falling outside the ± 5% error bands, indicating predictive instability. In contrast, CGEP2 (Figure 5b) shows superior robustness, as most points are clustered around the ideal diagonal and fit within ± 5% margins. Notably, although 3D printed concrete typically exhibits significant strength decay in the Z-direction (vertical loading) due to weak inter-layer interfaces, the median prediction error for CGEP2 in the Z-direction is only marginally higher than that in the X-direction. This indicates that CGEP2 accurately predicts the anisotropic behavior and accounts for the loss of strength due to inter-layer defects.

Display Full Size

Figure 5. Regression plots of predicted versus experimental values: (a) CGEP1 model; (b) CGEP2 model. CGEP: chaotic gene expression programming.

The observed performance disparity stems fundamentally from the physical applicability of the employed function sets. CGEP1 attempts to approximate highly non-linear material behaviors by relying solely on rational fractions constrained by basic arithmetic operations, which inevitably results in predictive dispersion. Conversely, incorporating non-linear operators into the function set of CGEP2 endows the model with superior mathematical expressiveness. Specifically, the Sqrt function facilitates the capture of non-linear characteristics inherent in fiber reinforcement, such as marginal effects, while the Ln function aligns closely with the kinetics of concrete strength evolution over time. CGEP2 effectively extracts latent physical laws from the data through the integration of these critical operators and optimized evolutionary parameters (e.g., higher mutation rates), thereby achieving exceptional predictive accuracy on the testing set.

Next, this study designed comparative experiments with two widely used models, multiple linear regression (MLR) and RF, to further validate the performance advantages of CGEP2. MLR represents traditional linear modeling approaches, while RF represents ensemble learning methods. All models were trained and evaluated on identical datasets to ensure fair comparison.

The performance metrics of the three models are summarized in Table 5. MLR achieved a testing R² of only 0.647 with an RMSE of 0.623 MPa, indicating that linear assumptions fail to capture the nonlinear behavior of fiber reinforcement and hydration evolution. RF, through an ensemble of 200 decision trees with optimized hyperparameters, achieved higher accuracy (testing R² = 0.883, RMSE = 0.359 MPa). However, its training R² (0.967) significantly exceeded its testing R², suggesting overfitting tendencies. More critically, the black-box nature of RF prevents engineers from obtaining explicit formulas for practical calculations.

In contrast, CGEP2 achieved testing R² = 0.914, RMSE = 0.312 MPa, and MAE = 0.244 MPa, outperforming RF by 3.5% in R² and 13.1% in RMSE reduction, while providing a fully interpretable analytical formula. Additionally, the minimal gap between training R² and testing R² demonstrates superior generalization capability without the overfitting issues observed in RF.

In summary, RF achieves high accuracy but sacrifices interpretability and exhibits overfitting, while MLR lacks predictive accuracy. In contrast, CGEP2 successfully balances both requirements, delivering superior performance with complete mathematical transparency.

3.1.2 Error distribution analysis

The absolute error distribution for the testing set (Figure 6) is utilized to scrutinize predictive stability at the individual sample level, elucidating point-wise error characteristics across the entire data domain.

Display Full Size

Figure 6. Distribution of absolute errors for CGEP predictions on the testing set. CGEP: chaotic gene expression programming.

The error trajectory for CGEP1 manifests pronounced instability (Figure 6), characterized by high-frequency, large-amplitude oscillations, where absolute error peaks for multiple samples surpass 0.8 MPa. Such erratic fluctuations suggest hypersensitivity to specific sample features and limited robustness. In contrast, CGEP2 exhibits exceptional stationarity, with errors consistently maintained within 0.3 MPa. The significantly reduced magnitude and frequency of errors for CGEP2, relative to CGEP1, indicate that this model not only excels in aggregate performance but also maintains superior consistency in local, individual sample predictions.

Overall, multi-dimensional evaluations based on R², RMSE, MAE, and error distribution confirm that CGEP2 significantly outperforms CGEP1 in predictive accuracy, generalization capability, and stability. CGEP2 effectively transcends the limitations of basic algebraic operations by leveraging an optimized gene topology and an extended non-linear function set, precisely capturing intricate non-linear coupling relationships between input variables and STS. Consequently, CGEP2 serves as the optimal model for subsequent sensitivity analysis and physical interpretation.

3.2 Sensitivity analysis

This section deconstructs the internal decision-making logic of the CGEP2 model to quantify the contribution weights of input parameters to STS. It elucidates the intrinsic influence of material and processing parameters on STS from dual dimensions of linear correlation and non-linear contribution through the evaluation of relevance factors and PFI.

3.2.1 Relevance factor analysis

Figure 7 shows the relevance of input variables to STS based on the training dataset. The magnitude of r signifies the strength of influence of each variable on STS, with positive and negative values denoting corresponding correlations.

Display Full Size

Figure 7. Relevance factor between input variables and predicted STS. CGEP: chaotic gene expression programming; FV: fiber volume; LH: layer height; LD: load direction; RCA-r: recycled aggregate replacement ratio.

The w/b ratio displays the highest negative correlation (r = -0.45), consistent with micromechanics theory, where higher water content increases capillary porosity and reduces matrix density. RCA-r also shows a significant negative correlation (r = -0.38), confirming the detrimental effect of attached porous mortar and weak interfacial transition zones. The printing parameters LD and LH exhibit distinct negative correlations (r = -0.25 and -0.21, respectively). The former captures the anisotropy where cracks propagate along weak inter-layer interfaces under vertical loading, while the latter indicates that increased layer thickness impairs bonding by reducing relative contact area or by enlarging pore defects. Conversely, the FV presents a positive correlation (r = 0.18), validating the bridging mechanism that effectively facilitates stress transfer and inhibits the transition from micro- to macro-cracks.

3.2.2 Permutation feature importance

A PFI analysis was performed to identify non-linear relationships of individual variables, in addition to the linear relationships defined by the relevance factors. Table 6 shows both PFI score and rank of importance for each input variable.

The PFI results align closely with the relevance factor analysis, confirming w/b and RCA-r as the paramount determinants of STS. This consensus indicates that matrix densification dictated by w/b and aggregate quality determined by the RCA replacement ratio constitute the fundamental intrinsic factors controlling tensile strength. Notably, the LD ranks third, surpassing both LH and FV. This ranking quantitatively underscores that the structural anisotropy of 3D printed concrete often exerts a governing influence on mechanical performance, outweighing the impact of specific compositional adjustments. Consequently, the presence of weak inter-layer interfaces represents a critical bottleneck limiting STS. Although LH and FV rank lower, their substantial PFI values suggest that rational regulation of layer thickness and fiber incorporation remain effective strategies for enhancing performance, provided that primary mix parameters are optimized.

In summary, the PFI analysis validates model robustness and offers explicit engineering guidance. For the mix design of 3D printed recycled concrete, priority should be assigned to controlling the w/b ratio and RCA content. Furthermore, rigorous alignment between printing paths and loading directions is essential to mitigate strength reduction induced by anisotropy.

3.2.3 Effect of 3D printing parameters and loading direction

This section presents a decoupling analysis of layer height and loading direction to elucidate the specific impacts of 3D printing parameters on mechanical performance.

The sensitivity analysis indicates a negative correlation between LH and STS (r = -0.21), reflecting the defect amplification effect of layer height. STS paradoxically declines as LH increases, despite the significant negative correlation between LH and the w/b ratios in the training dataset. In this dataset, higher layers typically correspond to lower w/b ratios (stiffer inks). This counter-intuitive trend reveals that physical defects induced by layer geometry override the benefits of matrix densification. Specifically, increased layer height enlarges the critical defect size and compromises inter-layer compaction. Consequently, these geometric defects negate the strength enhancement theoretically expected from lower w/b ratios.

The dependency of STS on loading direction represents the intrinsic distinction between 3D printed and cast concrete. The prediction error distribution for 3DPRAC samples across varying loading directions is depicted in Figure 8 to quantify the impact of this anisotropy on predictive accuracy. The statistical characteristics of the errors reveal distinct mechanisms governing prediction accuracy:

Display Full Size

Figure 8. Distribution of prediction errors for 3DPRAC samples under different loading directions. 3DPRAC: 3D printed recycled aggregate concrete.

(1) The median prediction error in the Z-direction (0.21 MPa) marginally exceeds that in the X-direction (0.19 MPa) and Y-direction (0.18 MPa). The high variability of interfacial bonding quality, influenced by factors such as printing time intervals and surface moisture evaporation, maximizes strength dispersion when tensile stress acts perpendicularly to the inter-layer interfaces. This consequently increases prediction difficulty.

(2) Although the X-direction (parallel to the printing path) exhibits a lower median error, its error distribution range (whiskers) is notably extensive, with a maximum error reaching 0.45 MPa. This is likely attributed to process-induced heterogeneity along the printing path. Nozzle start-stops, velocity fluctuations, and material accumulation at corners introduce additional microscopic defects in the X-direction, causing strength deviations in certain samples.

(3) In contrast, the Y-direction displays the most compact interquartile range and the narrowest error span. This indicates that the material exhibits the most stable mechanical response under lateral loading modes, resulting in the highest predictive confidence for the model.

In summary, the CGEP2 model not only successfully captures the strength anisotropy stemming from weak inter-layer interfaces (characterized by lower strength and higher unpredictability in the Z-direction) but also sensitively reveals process instability along the printing path (characterized by significant fluctuations in the X-direction).

3.3 Interpretable prediction formula

This study established an explicit formula applicable to both 3D printed and conventional recycled concrete, grounded in the evolutionary principles of GEP. The formula is constructed through the linear superposition of four Sub-ETs, as shown in Figure 9. The formula is expressed as:

Display Full Size

Figure 9. Sub-ETs structure of CGEP2 model. (Yellow nodes represent operators; pink nodes represent constants: c1 = 2.66, c2 = -1.19, c3 = 8.96, c4 = 1.02, c5 = 0.63, c6 = 1.59; green nodes represent input variables; arrows indicate the data flow from leaf nodes to root nodes). CGEP: chaotic gene expression programming; STS: splitting tensile strength.

(7)$STS=\sqrt{2.66\cdot x_{FV}}-1.19\cdot x_{LD}+\frac{8.96-x_{w/b}\cdot \sqrt{x_{RCA-r}}}{1.02\cdot x_{LH}+\frac{x_w}{b}+0.63}+1.59\cdot \ln \left(\mathrm{Age}\right)$

where x_w/b represents the water-to-binder ratio, x_RCA-r the recycled concrete aggregate replacement ratio, x_FV the fiber volume fraction, x_LH the layer height, x_LD the loading direction, and Age the curing period.

Each term in Eq. (7) corresponds to an independently evolving sub-tree structure that captures a specific physical mechanism, with the linear combination of all sub-trees constituting the complete predictive model.

The initial term $\sqrt{2.66\cdot x_{FV}}$ (Sub-ET1 in Figure 9) captures the non-linear toughening effect of fibers on the matrix. The square root operator in this sub-tree structure simulates the diminishing marginal returns of fiber reinforcement. At low dosages, fibers disperse uniformly with high bridging efficiency; at high dosages, fiber agglomeration reduces the marginal reinforcement effect.

The linear term -1.19·x_LD (Sub-ET2 in Figure 9) serves as an anisotropic penalty term in the model to quantify the loading direction effect of 3DCP. This term is a binary variable that equals 0 for parallel loading and 1 for perpendicular loading, imposing a strength penalty of 1.19 MPa. Perpendicular loading applies principal tensile stress directly to the weak interlayer interfaces, activating microscopic defects at the interfaces, and triggering strength degradation. The magnitude of this effect is primarily determined by the bonding quality of interlayer interfaces and exhibits minimal interaction with other parameters such as w/b ratio and RCA-r, thus manifesting as an independent linear penalty term.

Subsequently, the fractional term $\frac{8.96-x_{w/b}\cdot \sqrt{x_{RCA-r}}}{1.02\cdot x_{LH}+x_{\frac{w}{b}}+0.63}$ (Sub-ET3 in Figure 9) is the most complex in structure among the subtrees and represents the coupling term between material composition and printing parameters. The numerator term, representing the interaction between x_w/b and x_RCA-r, quantifies the strength deterioration induced by recycled aggregates. This product term reflects the significant synergistic degradation effect between the w/b ratio and the recycled concrete aggregate replacement ratio. The simultaneous increase of two parameters leads to a nonlinear amplification of the negative impact on STS. Quantitative analysis reveals that when w/b = 0.60 and RCA-r = 100%, the strength reduction reaches 20.7%, which far exceeds the 7.3% strength reduction caused by the increase of RCA-r alone, confirming the existence of a nonlinear synergistic degradation mechanism. The denominator term captures the size effect of layer height, where the positive coefficient (1.02) of x_LH indicates that as the layer height increases, the sizes of voids and defects at the interlayer interface increase, leading to a decrease in interlayer compaction and an increase in porosity. This continuous and gradual effect aligns with the scale effect law of geometric defects.

Finally, the temporal evolution of strength is modeled by the hydration age correction term +1.59·ln(Age) (Sub-ET4 in Figure 9). The logarithmic form aligns with the reaction kinetics of concrete hydration, where strength develops rapidly in the early stages but tends to plateau as the age extends and the hydration rate decelerates. This term provides a baseline correction to reflect the time-evolving maturation process of the cementitious matrix.

The interpretable STS prediction model was developed based on the dataset with specific parameter ranges. The validated parameter ranges based on the dataset are shown in Table 7.

3.4 Limitations and future perspectives

The CGEP2 model is accurate and interpretable within its validated range. The model is demonstrably applicable in predicting the compressive strength of 3D printed concrete, fiber-reinforced concrete, recycled aggregate concrete, and conventional concrete. This enables rapid strength prediction at the preliminary design stage, making it useful for mix proportioning and performance evaluation in engineering practice. However, further assessment is required to refine or constrain its applicability.

The current dataset concerns intrinsic material attributes, where traditional concrete is isotropic while 3D printed concrete exhibits anisotropy. The LD variable is tied to the 3DPRAC database, creating an inherent coupling between anisotropy and the printing process. Future research incorporating diverse printing paths (e.g., concentric and cross-hatching) would facilitate a more granular decoupling of geometric path effects from material composition.

The derived explicit formula currently isolates layer height as the primary geometric constraint. However, mechanical performance is fundamentally governed by the interplay between rheology and hardening kinetics, requiring consideration of additional parameters such as printing speed, nozzle diameter, and time interval. Future studies should incorporate these process variables to enhance model generalization across diverse printing scenarios.

The model construction in this study relies on 110 experimental data sets. While this sample size proves sufficient for the GEP algorithm to effectively capture key physical laws, incorporating more data would further strengthen the generalization capability of the model across diverse material systems and printing conditions. This dataset expansion, combined with emerging techniques such as transfer learning aligned with the Materials Genome Initiative, offers promising pathways to enhance predictive performance. Fine-tuning pre-trained models with incrementally expanded datasets may help transcend current data scale limitations while maintaining computational efficiency.

4. Conclusions

In this study, a unified GEP framework to predict concrete splitting tensile strengths was created. This model was developed based on 110 datasets, including RAC, fiber-reinforced RAC, and 3DPRAC. The main conclusions are as follows:

(1) The developed CGEP2 model demonstrated exceptional robustness, achieving an R² of 0.914 and an RMSE of 0.312 MPa. The model successfully integrated distinct material systems by decoupling material types into continuous physical descriptors (e.g., RCA-r and FV), capturing nonlinear relationships between microscopic components and macroscopic process parameters within a unified computational framework.

(2) Sensitivity analysis established that LD and LH are the key drivers for mechanical performance with respect to 3D printed concrete, with their influence outweighing material constituents such as fiber volume. The model revealed that loading perpendicular to the layers (Z-direction) induces significant strength degradation, confirming the weak interlayer interface mechanism and suggesting the need for contact optimization of printing paths.

(3) An interpretable STS prediction equation was constructed to describe the saturation of fiber content, layer height defect, and anisotropy penalties as mathematical terms. The explicit equations provide a rapid, low-cost tool for mix design and strength prediction, replacing black-box software with a tool that engineers can use directly.

Despite the high fidelity of the developed framework, future work should focus on three extensions. First, integrating diverse printing paths (e.g., concentric and cross-hatching) is necessary to decouple geometric path effects from material composition. Second, incorporating other high-dimensional process variables (e.g., printing speed and interval time) increases the chance for model generalization to more complex printing scenarios. Finally, transfer learning techniques offer a pathway to overcome current data limitations and further improve predictive robustness.

Acknowledgements

The authors declare that AI tools were used solely for language polishing during the manuscript preparation process. All research content, including study design, data analysis, interpretations, figures, and tables, is original and was not generated using AI tools.

Authors contribution

Wang T: Conceptualization, methodology, formal analysis, writing-original draft.

Liu C: Conceptualization, methodology, writing-review & editing.

Zhu C: Formal analysis, writing-review & editing.

Zhu W: Writing-review & editing.

Liu H: Investigation, writing-review & editing.

Conflicts of interest

Huawei Liu serves as a Youth Editorial Board Member of the Journal of Building Design and Environment. The other authors declare no conflicts of interest.

Ethical approval

Not applicable.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Funding

This study was supported by the National Natural Science Foundation of China (52178251 and 62276207), the National Natural Science Foundation of China Youth Fund (52408289 and 52508304), the Technology Innovation Guidance Program of Shaanxi Province (2023GXLH-049), Shaanxi Provincial Natural Science Basic Research Program (2025JC-YBMS-550), Qinchuangyuan’s “Scientist and Engineer” Team Building of Shaanxi Province (2023KXJ-242), and the Special Research Program for Local Service of Shaanxi Province (23JC047)

Copyright

© The Author(s) 2026.