Beams as structural components are crucial in many structural systems. The prediction of their deflection is essential since excessive values can lead to the structural system losing its operational serviceability (Serviceability Limit State—SLS). On the other hand, composite materials are increasingly used in structural engineering due to their enhanced stiffness combined with reduced weight. Several shear deformation theories have been developed so far to evaluate the response of thin, moderately thick, or deep beams. They fall into three main categories: the Euler-Bernoulli beam theory (or Classical Beam Theory—CBT), the Timoshenko beam theory (or First Order Beam Theory—FOBT) and the Higher-Order Beam Theories (HOBTs). CBT is applicable for thin beams with no shear effect. In the FOBT, a constant state of transverse shear strain is assumed that does not satisfy the zero shear stress condition at the top and bottom edges of the beam and thus requires a shear correction factor to compensate for this error (see, e.g., Wang et al., 2000; Eisenberger, 2003; Civalek and Kiracioglu, 2010; Lin and Zhang, 2011; Endo, 2016). In general, the HOBTs adopt a specific function (parabolic, trigonometric, exponential, or hyperbolic) to more accurately represent the shear stress distribution along the beam’s thickness and do not require the shear correction factor (see e.g., Reddy, 1984; Heyliger and Reddy, 1988; Khdeir and Reddy, 1997; Murthy et al., 2005; Vo and Thai, 2012; Pawar et al., 2015; Nguyen et al., 2017; Srinivasan et al., 2019). The literature contains a plethora of publications on the subject, and the interested reader is referred to the excellent review paper of Liew et al. (2019). In this investigation, a refined higher-order beam theory is utilized for the analysis of laminated composite beams based on Reddy-Bickford’s third-order beam theory (Wang et al., 2000) which was derived independently by Bickford (1982) and Reddy (1984).

Utilizing higher-order beam theories for more accurate analyses entails a significant increase in complexity as compared to low-order theories, as the latter are mathematically simpler and more widely used. The main motivation of this work is to bridge this gap and provide a simple computational tool to allow for the fast design of beams while keeping the best of both worlds, i.e., the more accurate results of a refined high-order theory and the ease of application of the low-order theories. In order to achieve that, the geometric and material variables are discretized within fairly wide, yet reasonable ranges. After applying the high-order analyses, the results are collected, tabulated, and used as input for multiple machine learning algorithms, i.e., regression models. These models provide a fast and easy-to-use computational tool that can be used for preliminary design and optimization. Regression analysis also yields important insights regarding the performance of each model, the effect of boundary conditions, and the relative importance of each input variable for the problem at hand.

The rest of the paper is organized as follows. A theoretical formulation of the problem is carried out and explained in detail next, followed by a summary of the regression methods utilized in this work. The numerical results are presented next, along with their discussion. Finally, the conclusions drawn based on the findings of this work are presented.

Consider an elastic symmetric cross-ply laminated rectangular beam ( b × h ) of length l , with x being the axial coordinate and z being the coordinate along the thickness of the beam. The fibers of each ply are aligned at an angle θ with respect to the x axis (see Figure 1).

The beam is subjected to a transverse distributed loading p z , respectively. Based on the higher-order theory for laminated composite plates introduced by Reddy (1984), the displacement field of an arbitrary point on the beam cross-section is given by u 1 ( x , z ) = z [ ψ ( x ) − 4 3 z 2 h 2 ( ψ ( x ) + ∂ w ( x ) ∂ x ) ] (1) u 2 ( x , z ) = 0 (2) u 3 ( x , z ) = w ( x ) (3) where w ( x ) is the transverse displacement of the midplane ( z = 0 ); ψ ( x ) is the rotation of a normal to the midplane, and x , z are the axial and thickness coordinates of the beam.

Splitting the transverse displacement w ( x ) into a bending w b ( x ) and a shear w s ( x ) component, i.e., Vo and Thai (2012). w ( x ) = w b ( x ) + w s ( x ) (4) and introducing the transformation ∂ w s ( x ) ∂ x = ψ ( x ) + ∂ w ( x ) ∂ x   or   ψ ( x ) = − [ ∂ w ( x ) ∂ x − ∂ w s ( x ) ∂ x ] = − ∂ w b ( x ) ∂ x (5)

Equations 1–3 can be rewritten in the following form u 1 ( x , z ) = − z ∂ w b ( x ) ∂ x − f ( z ) ∂ w s ( x ) ∂ x (6) u 2 ( x , z ) = 0 (7) u 3 ( x , z ) = w b ( x ) + w s ( x ) (8) where f ( z ) = 4 3 z 3 h 2 . The displacement field given above yields the following nonzero components of the strain tensor ε x = − z ∂ 2 w b ∂ x 2 − f ( z ) ∂ 2 w s ∂ x 2 (9) γ x z = [ 1 − d f ( z ) d z ] ∂ w s ∂ x = g ( z ) ∂ w s ∂ x (10) where g ( z ) = ( 1 − 4 z 2 h 2 ) , and for reasons of brevity w b = w b ( x ) and w s = w s ( x ) .

Substituting Eqs 9, 10 into the stress-strain relations for the kth lamina in the lamina coordinate we obtain (Khdeir and Reddy, 1997) σ x ( k ) = Q ¯ 11 ( k ) ε x (11) τ x z ( k ) = Q ¯ 55 ( k ) γ x z (12) with Q ¯ 11 ( k ) , Q ¯ 55 ( k ) being the well-known transformed elastic stiffnesses Q ¯ 11 ( k ) = Q 11 ( k ) ⁡ cos 4 θ k + 2 ( Q 12 ( k ) + 2 Q 66 ( k ) ) sin 2 θ k ⁡ cos 2 θ k + Q 22 ( k ) ⁡ sin 4 θ k (13) Q ¯ 55 ( k ) = Q 44 ( k ) ⁡ sin 2 θ k + Q 55 ( k ) ⁡ cos 2 θ k (14) and Q 11 ( k ) , Q 12 ( k ) , Q 22 ( k ) , Q 44 ( k ) , and Q 55 ( k ) are Q 11 ( k ) = E 1 ( k ) 1 − ν 12 ( k ) ν 21 ( k ) ,   Q 12 ( k ) = E 2 ( k ) ν 12 ( k ) 1 − ν 12 ( k ) ν 21 ( k ) ,   Q 22 ( k ) = E 2 ( k ) 1 − ν 12 ( k ) ν 21 ( k ) (15) Q 44 ( k ) = G 23 ( k ) ,   Q 55 ( k ) = G 13 ( k ) ,   Q 66 ( k ) = G 12 ( k ) (16) while θ k is the angle between the principal material axis and the coordinate x axis.

Applying the Principle of Virtual Work ∫ 0 l ∫ A [ σ x ( k ) δ ε x + τ x z ( k ) δ γ x z ] d A d x − ∫ 0 l [ p z δ ( w b + w s ) ] d x = 0 (17) and substituting Eqs 9, 10 yields ∫ 0 l ∫ A { σ x ( k ) [ − z ∂ 2 δ w b ∂ x 2 − f ( z ) ∂ 2 δ w s ∂ x 2 ] + τ x z ( k ) g ( z ) ∂ δ w s ∂ x } d A d x − ∫ 0 l [ p z ( δ w b + δ w s ) ] d x = 0 (18)

Introducing now the following stress resultants M b = ∫ A z σ x ( k ) d A ,   M s = ∫ A f ( z ) σ x ( k ) d A ,   Q = ∫ A g ( z ) τ x z ( k ) d A (19) Eq. 18 become ∫ 0 l ( − M b ∂ 2 δ w b ∂ x 2 − M s ∂ 2 δ w s ∂ x 2 + Q ∂ δ w s ∂ x ) d x − ∫ 0 l [ p z ( δ w b + δ w s ) ] d x = 0 (20)

Integrating the appropriate terms in the above equation and collecting the coefficients of δ w b , and δ w s we obtain the following governing equations ∂ 2 M b ∂ x 2 = − p z (21) ∂ 2 M s ∂ x 2 + ∂ Q ∂ x = − p z (22) together with the following associated boundary conditions of the form: specify w b   or   Q b ≡ ∂ M b ∂ x (23) w s   or   Q s ≡ ∂ M s ∂ x + Q (24) ∂ w b ∂ x   or   M b (25) ∂ w s ∂ x   or   M s (26)

Substituting Eqs 11, 12 into Eq. 19 and using Eqs 9, 10 yields the stress resultants in terms of the displacements as M b = − D 11 ∂ 2 w b ∂ x 2 − F 11 ∂ 2 w s ∂ x 2 (27) M s = − F 11 ∂ 2 w b ∂ x 2 − H 11 ∂ 2 w s ∂ x 2 ,   Q = A 55 ∂ w s ∂ x (28) where D 11 = b ∫ − h / 2 h / 2 Q ¯ 11 ( k ) z 2 d z ,   F 11 = b ∫ − h / 2 h / 2 Q ¯ 11 ( k ) z f ( z ) d z ,   H 11 = b ∫ − h / 2 h / 2 Q ¯ 11 ( k ) f 2 ( z ) d z (29) A 55 = b ∫ − h / 2 h / 2 Q ¯ 55 ( k ) g 2 ( z ) d z (30)

Finally, after the substitution of the stress resultants, Eqs 27, 28 into Eqs 21, 22, we arrive at the equilibrium equations in terms of the displacements − D 11 ∂ 4 w b ∂ x 4 − F 11 ∂ 4 w s ∂ x 4 = − p z (31) − F 11 ∂ 4 w b ∂ x 4 − H 11 ∂ 4 w s ∂ x 4 + A 55 ∂ 2 w s ∂ x 2 = − p z (32) which together with the pertinent boundary conditions (23)–(26) constitute the boundary value problem solved using the Analog Equation Method (AEM), a robust numerical method based on an integral equation technique (Katsikadelis and Tsiatas, 2003; Tsiatas et al., 2018).

In this work, several linear and nonlinear regression models are comparatively examined. Linear regression is a linear model that assumes a linear relationship between the input variables and the output variable, and the predicted value can be calculated from a linear combination of the input variables (Narula and Wellington, 1982). The distance from each data point to the predicted values is calculated and sum all these squared errors together. This quantity is minimized by the ordinary least squares method to estimate the optimal values for the coefficients of each independent variable.

There are extensions of the linear model called regularization methods. These methods seek to both minimize the sum of the squared error of the model on the training set but also to reduce the complexity of the model. Two popular regularization methods for linear regression are the Lasso Regression (Zou et al., 2007) where Ordinary Least Squares is modified to also minimize the absolute sum of the coefficients (L1 regularization), and the Ridge Regression (Hoerl et al., 1985) where Ordinary Least Squares is modified to also minimize the squared absolute sum of the coefficients (L2 regularization). A Bayesian view of ridge regression is obtained by noting that the minimizer can be considered as the posterior mean of a model (Tipping, 2001). The elastic net (Friedman et al., 2010) is a regularized regression method that linearly combines the L1 and L2 penalties of the lasso and ridge methods. Huber’s criterion is a hybrid of squared error for relatively small errors and absolute error for relatively large ones. Lambert-Lacroix and Zwald (2011) proposed Huber regressor to combine Huber’s criterion with concomitant scale and Lasso.

An L1 penalty minimizes the size of all coefficients and allows any coefficient to go to the value of zero, acting as a type of feature selection method since removes input features from the model. Least Angle Regression (Efron et al., 2004) is a forward stepwise version of feature selection for regression that can be adapted for the Lasso not to require a hyperparameter that controls the weighting of the penalty in the loss function since the weighting is discovered automatically by Least Angle Regression method via cross-validation. LassoLars is a lasso model implemented using the Least Angle Regression algorithm, where unlike the implementation based on coordinate descent, this yields the exact solution, which is piecewise linear as a function of the norm of its coefficients.

Orthogonal matching pursuit (Pati et al., 1993) tries to find the solution for the L0-norm minimization problem, while Least Angle Regression solves the L1-norm minimization problem. Although these methods solve different minimization problems, they both depend on a greedy framework. They start from an all-zero solution, and then iteratively construct a sparse solution based on the correlation between features of the training set and the output variable. They converge to the final solution when the norm approaches zero.

K Neighbors Regressor (KNN) algorithm uses feature similarity to predict the values of new instances (Altman, 1992). The distance between the new instance and each training instance is calculated, the closest k instances are selected based on the preferred distance and finally, the prediction for the new instance is the average value of the dependent variable of these k instances.

Unlike linear regression, Classification and Regression Tree (CART) does not create a prediction equation, but data are partitioned into subsets at each node according to homogeneous values of the dependent variable and a decision tree is built to be used for making predictions about new instances (Breiman et al., 1984). We can enlarge the tree until always gives the correct value in the training set. However, this tree would overfit the data and not generalize well to new data. The correct policy is to use some combination of a minimum number of instances in a tree node and maximum depth of tree to avoid overfitting.

The basic idea of Boosting is to combine several weak learners into a stronger one. AdaBoost (Freund and Schapire, 1997) fits a regression tree on the training set and then retrains a new regression tree on the same dataset but the weights of each instance are adjusted according to the error of the previous tree predictions. In this way, subsequent regressors focus more on difficult instances.

Random Forests algorithm (Breiman, 2001) builds several trees with the CART algorithm using for each tree a bootstrap replica of the training set with a modification. At each test node, the optimal split is derived by searching a random subset of size K of candidate features without replacement from the full feature set.

Like Random Forests, Gradient Boosting (Friedman, 2001) is an ensemble of trees, however, there are two main differences. Firstly, the Random forests algorithm builds each tree independently while Gradient Boosting builds one tree at a time since it works in a forward stage-wise manner, introducing a weak learner to improve the shortcomings of existing weak learners. Secondly, Random Forests combine results at the end (by averaging the result of each tree) while Gradient Boosting combines results during the process.

LightGBM (Ke et al., 2017) extends the gradient boosting algorithm by adding automatic feature selection and focusing on instances with larger gradients to speed up training and sometimes even improve predictive performance.

The Extra-Trees algorithm (Geurts et al., 2006) creates an ensemble of unpruned regression trees according to the well-known top-down procedure of the regression trees. The main differences concerning other tree-based ensemble methods are that the Extra-Trees algorithm splits nodes by choosing fully at random cut-points and that uses the whole learning set (instead of a bootstrap replica) to grow the trees.

Passive-Aggressive regressor (Crammer et al., 2006) is generally used for large-scale learning since it is an online learning algorithm. In online learning, the input data come sequentially, and the learning model is updated step-by-step, as opposed to batch learning, where the entire dataset is used at once.

The scope of the current study is to exploit predictive models for the maximum deflection max w of a symmetric cross-ply ( θ 1 ° / θ 2 ° / θ 3 ° ) rectangular beam for various span-to-depth ratios and boundary conditions subjected to a uniformly distributed load p z . All laminates are of equal thickness and made of the same orthotropic material. The main parameters that influence the response of the composite beams are the moduli of elasticity E 1 , E 2 , the span-to-depth L = l / h and the ply angles θ 1 ° / θ 2 ° / θ 3 ° . The range of values of the parameters together with the material properties are given as: E 1 , E 2 = { 1,2,3 , ...15 } , G 12 = G 13 = 0.5 E 2 , G 23 = 0.2 E 2 , ν 12 = 0.25 , L = l / h = { 1,2,3,4,5 } , θ 1 , θ 2 = { 0 , π / 8 , π / 4 ,   3 π / 4 , π / 2 } , and θ 1 = θ 3 (due to symmetry). For the given range of the parameters, Eqs 31, 32 are solved numerically producing a comprehensive database for each one of the examined boundary conditions presented in Table 1. This dataset contains 15 × 15 × 5 × 5 × 5 = 28125 values of max w which are used in the regression analysis.

A plethora of regression algorithms, presented in the previous section, were employed for building corresponding predictive models of the max w using pyCaret (Ali, 2020), which is an open-source software machine learning library. A 5-fold cross-validation resampling procedure was used for evaluating the performance of the predictive models. The dataset was randomly divided into five folds of equal size and each fold was used for evaluating the performance of the model trained on the rest folds, whereas the final measure was the average value of the computed evaluation metrics on each test fold. Evaluation metrics are a measure of how well a model performs. The most popularly used evaluation metrics for regression problems are the mean absolute error (MAE), the mean absolute percentage error (MAPE), the mean square error (MSE), the root mean square error (RMSE), the root mean squared log error (RMSLE) and the coefficient of determination. The lower the value of these metrics the better the model. The perfect value of metrics is 0, indicating that the prediction model is perfect. To quantify the accuracy of the examined algorithms, the following evaluation metrics are used herein:

Mean absolute error (MAE) MAE = 1 n ∑ i = 1 n | y ′ i − y i | (33)

Mean absolute percentage error (MAPE) MAPE = 1 n ∑ i = 1 n ( | y i ′ − y i | ) / | y i ′ | (34)

Mean square error (MSE) MSE = 1 n ∑ i = 1 n ( y i ′ − y i ) 2 (35)

Root mean square error ( RMSE ) RMSE = 1 n ∑ i = 1 n ( y i ′ − y i ) 2 (36)

Root mean squared log error ( RMSLE ) RMSLE = 1 n ∑ i = 1 n ( log ( y i ′ + 1 ) − log ( y i + 1 ) ) 2 (37)

Coefficient of determination ( R 2 ) R 2 = S S r e g S S t o t = ∑ j ( y i ′ − y ¯ ) 2 ∑ i ( y i − y ¯ ) 2 , (38) where y ′ i refers to predicted values, and y i refers to true values. S S r e g is the regression sum of squares (i.e., explained sum of squares), and S S t o t is the total sum of squares, which is proportional to the variance of the data. The coefficient of determination ( R 2 ) is the square of the correlation between the actual and predicted variable and ranges from 0 to 1 . A zero value indicates that the model cannot explain any of the predicted variables. A value of 1 indicates that the regression model explains perfectly the predicted variable.

Apart from the evaluation metrics of the machine learning algorithms, two other useful tools are presented for the predictive analysis of the max w . First, the feature importance is a technique for assigning scores to input features that indicate the relative importance of each feature for the prediction. The scores can highlight which features are most relevant to the target and the opposite, i.e., which features are the least relevant. Most importance scores are calculated using the most accurate predictive model that has been fit on our data (Louppe et al., 2013). Second, the correlation matrix heatmap illustrates the correlation dependence between the variables of the database. That is, each square of the matrix represents the correlation between the attributes paired on the two axes. A value of + 1 (or − 1 ) indicates a perfect correlation between two variables, with + 1 indicating a positive correlation and − 1 a negative (inverse) correlation; a value in the range from 0.6 to 1 (or from − 0.6 to − 1 ) indicates a strong correlation; a value between 0.4 and 0.6 (or between − 0.4 and − 0.6 ) indicates a moderate correlation; a value in the range from 0 to 0.4 (or from 0 to − 0.4 ) indicates a weak correlation.

In this paper, several machine learning regression models were employed for the prediction of the deflection of symmetric laminated composite beams subjected to a uniformly distributed load. Training, validation, and testing of the models require large amounts of data that cannot be provided by the scarce experiments. Instead, ample amounts of data are generated numerically using a refined higher-order beam theory for various span-to-depth ratios and boundary conditions, by appropriate discretization of all pertinent geometric and material properties.

The main conclusion that can be drawn from this investigation are as follows:

• Regarding the regression models, the Extra-Trees algorithm is, without doubt, the best performer for all cases of boundary conditions, followed by the Random Forest Regressor, the Decision Tree Regressor, the Light Gradient Boosting Machine, and the K Neighbors Regressor.