Hybrid analytical–machine-learning surrogate modeling for buckling prediction and design of functionally graded porous nanobeams

Functionally graded porous (FGP) structures combine low weight with high strength, energy absorption capacity, and favorable thermal and acoustic insulation. Their mechanical and functional properties can be tailored by controlling the size, density, gradation, and spatial distribution of pores [1], [2], [3]. These characteristics have supported the use of FGP structures in aerospace, automotive, civil, and biomedical engineering. Their increasing engineering relevance has also led to sustained interest in modeling their mechanical behavior under different loading and operating conditions. In particular, the buckling response of FGP beams has been examined using a range of structural theories and solution procedures. Magnucki and Stasiewicz [4] studied the buckling of porous beams under compressive loading using a broken-line hypothesis and the total potential energy principle. Magnucka-Blandzi and Magnucki [5] investigated the bending and buckling behavior of simply supported sandwich beams with metal foam cores and considered the associated design problem through two different approaches. Chen et al. [6] analyzed the bending and buckling of shear-deformable FGP beams using Timoshenko beam theory and the Ritz method, showing that symmetric porosity distributions produced higher buckling loads and smaller transverse deflections than asymmetric distributions. Chen et al. [7] subsequently examined the nonlinear vibration and post-buckling behavior of metal-foam FGP beams reinforced with graphene platelets. Kitipornchai et al. [8] investigated the free vibration and elastic buckling of multilayered functionally graded (FG) metal-foam beams reinforced with graphene platelets. Nguyen et al. [9] developed solutions for the bending, buckling, and free vibration of shear-deformable FGP beams using Lagrange’s principle and the Ritz method. Wu et al. [10] compared classical and shear-deformation beam theories for the buckling and free-vibration analysis of FGP beams through the generalized differential quadrature method.

Advances in nanotechnology have enabled the development of nanoelectromechanical systems (NEMS) with distinctive mechanical, electrical, and thermal characteristics. These systems are used in medical, optical, electronic, and energy-related applications and commonly contain nanoscale structural components such as beams, plates, and shells [11]. Nanobeams serve as one-dimensional load-carrying or sensing components in devices such as atomic-force-microscope cantilevers, nanoactuators, resonators, and nanosensors. The performance and accuracy of such devices therefore depend strongly on the mechanical response of their nanobeam components [12], [13], [14], [15], [16].

At the nanoscale, structural behavior may differ substantially from that predicted for macroscale components because conventional continuum theories do not account for size-dependent effects. Nonlocal strain-gradient theory, modified couple-stress theory, and nonlocal elasticity theory have consequently been developed or adopted to represent these effects. Among these formulations, Eringen’s nonlocal elasticity theory remains widely used in nanobeam analysis because it introduces the small-scale effect through a relatively concise constitutive framework [17], [18], [19], [20], [21]. The mechanical behavior of nanobeams has been investigated through several continuum formulations and numerical solution techniques [22], [23], [24]. Reddy [25] derived nonlocal governing equations for several beam theories and used a finite-element formulation to study the bending, buckling, and free vibration of nanobeams. Aydoğdu [26] developed a generalized nonlocal beam model for their static and dynamic analysis. Zhang et al. [24] proposed a two-parameter hybrid nonlocal beam model to represent the static and dynamic responses of micro- and nanobeams. Thai [27] and Thai and Vo [28] combined nonlocal elasticity with different shear-deformation theories to investigate nanobeam bending, buckling, and free vibration. Ebrahimi and Salari [29] examined the free vibration and thermal buckling of FG Timoshenko nanobeams using a Navier solution. Civalek and Demir [30] modeled microtubules as nonlocal Euler–Bernoulli beams embedded in an elastic medium and evaluated their buckling response using the finite-element method. Eltaher et al. [31] also employed a finite-element formulation to study the bending and vibration of FG Euler–Bernoulli nanobeams. Aria and Friswell [32] developed a size-dependent Timoshenko beam model within the nonlocal elasticity framework for the finite-element analysis of buckling and free vibration in FG nanobeams.

Reliable prediction of static buckling in FGP nanobeams is important for the design of stable NEMS components because their response is governed jointly by stiffness gradation, porosity distribution, beam geometry, and size-dependent behavior [8], [31]. Shafiei and Kazemi [33] investigated the buckling of two-dimensional tapered FGP nano- and microbeams by combining nonlocal elasticity with modified couple-stress theory. Mirjavadi et al. [34] used the generalized differential quadrature method to analyze the buckling and nonlinear vibration of FGP nanobeams. Aria et al. [35] studied the free vibration and buckling responses of FGP nanobeams subjected to thermal loading. More recently, Madasamy and Wang [36] examined the buckling behavior of bidirectional FGP nanobeams under moving loads in hygrothermal environments using a higher-order Haar wavelet method together with the generalized differential quadrature method. These studies have established useful analytical and computational foundations, but evaluating numerous combinations of geometric, material, and nonlocal parameters remains demanding when the analysis is extended to repeated prediction, uncertainty assessment, and design optimization.

Machine-learning (ML) methods provide an alternative means of constructing surrogate representations of structural responses across multidimensional parameter spaces [37]. Their use in the modeling of micro- and nanoscale structures has increased in recent years, particularly where repeated evaluations are required [38], [39], [40]. Tran et al. [41] combined an artificial neural network (ANN) with an optimization procedure to predict the buckling and vibration responses of FGP microplates. Tariq et al. [42] applied boosting-based ensemble models to investigate the effects of multiple parameters on the vibration response of porous nanobeams and used SHapley Additive exPlanations (SHAP) to assess the contribution of individual inputs. Al-Dahidi et al. [43] compared regression models, ANNs, and support vector machines for predicting the buckling response of FGP nanobeams under different boundary conditions. Cheng et al. [44] used an ANN to estimate the effects of eight parameters on the vibration behavior of FGP nanobeams. Although these studies demonstrate the predictive capacity of ML, most treat the learning model primarily as a direct input–output approximation. Less attention has been paid to the structured integration of an analytical mechanics formulation, residual learning, validation without information leakage, uncertainty quantification, interpretable parameter assessment, and multi-objective engineering design within a single framework.

The present study addresses this methodological gap by developing a hybrid analytical–ML surrogate framework for the buckling prediction and design of FGP nanobeams. The analytical component was formulated using exponential shear deformation theory, with three porosity distributions defined according to the Gibson–Ashby model. The exponential shear function satisfies the traction-free conditions at the upper and lower surfaces, removing the need for a shear correction factor. Eringen’s nonlocal elasticity theory was incorporated to represent small-scale effects, and the governing equations were derived from the minimum total potential energy principle. Navier’s method was then used to obtain the buckling loads of simply supported FGP nanobeams, and the formulation was verified against published benchmark results. Parametric analyses were conducted to determine the effects of the nonlocal parameter, slenderness ratio, porosity coefficient, and porosity distribution on the dimensionless buckling load (DBL). The resulting analytical dataset was subsequently used to construct a residual-learning surrogate in which a regularized mechanics-informed baseline represented the principal response trend, and nonlinear ML models captured the remaining discrepancy. The integrated framework further combined nested cross-validation (CV), conformal uncertainty quantification, SHAP-based interpretation, and Nondominated Sorting Genetic Algorithm II (NSGA-II) multi-objective optimization. This arrangement links analytical mechanics, predictive modeling, uncertainty assessment, model interpretation, and design exploration while preserving the physical trends of the underlying buckling formulation.

The FGP nanobeam was considered with length ($L$), thickness ($h$), and width ($b$), which were located on the x, y, and z axes, respectively. The beam was subjected to an axial load $N_0$ as shown in Figure 1a. The FGP nanobeam ’s material properties were assumed to either change or remain constant along the thickness direction, depending on the specific porosity distribution types (PDTs). Specifically, the functionally graded X-pattern (FGX) and functionally graded O-pattern (FGO) configurations denote a continuous decrease or increase in porosity from the mid-plane of the nanobeam to its top and bottom surfaces, respectively. Uniform distribution (UD) indicates that the porosity is distributed uniformly [10]. The effective Young’s and shear moduli can be defined as follows, respectively:

where, $v$, $E_{\max}$, and $G_{\max}$, represent the Possion’s ratio, the maximum values of the Young’s and shear moduli of the FGP nanobeam, respectively, as well as $\eta(z)$ and $\alpha_0$ represent the porosity distribution function and porosity coefficient, respectively, and are defined as follow:

where, $\eta_0$ is the PDT function for the UD, $E_{\min}$ and $G_{\min}$ are the minimum values of Young’s and shear moduli of the layer with the highest porosity, respectively.

According to the Gibson-Ashby micromechanical model, the relation between Young’s modulus and mass density is defined as \citep{45}:

Here, $\rho_{\max }$ is the mass density of non-porous material. The effective mass density of the FGP nanobeam can be formulated by substituting Eq. (4) into Eq. (1),

FGP nanobeam with distinct PDTs must have equal total mass, thus $\eta_0$ can be acquired:

The porosity coefficient for mass density $\alpha_\text{m}$ can be defined as:

where, $\rho_{\max}$ and $\rho_{\min }$ are the maximum and the minimum values of mass density, respectively.

Based on the present theory, axial displacement $u$ and transverse displacement $w$ of any point of the FGP nanobeam can be given as:

where, $u_0$ and $w_0$ are axial and transverse displacements of any point on the midplane of the FGP nanobeam, $\varphi$ is shear slope. $f(z)$ is the exponential shear-strain function which satisfies the traction-free boundary conditions on the nanobeam's outer surfaces \citep{46}. The strain components associated with the displacement field are:

where $\varepsilon_x$, $\gamma_{x z}$ and $g(z)$ are the axial strain, transverse shear strain and the first derivative of the shear-strain function, respectively.

The normal $\sigma_x$ and transverse $\tau_{x z}$ shear stresses at any point in a nanobeam can be obtained using Hooke’s Law \citep{3}:

According to the principle of minimum total potential energy, the first variation of the total potential energy of the system must be zero. Thus, governing equations of the FGP nanobeam can be derived as follows \citep{47}:

where, $\delta U$ and $\delta V$ are the first variation of the strain energy and the work done by external forces, respectively.

The first variation of the strain energy of the FGP nanobeam can be stated as:

where, $N$, $M^b$, $M^s$ and $Q$ are force and moment resultants which can be obtained by integrating stress components over the cross-sectional area as follows:

here, $A_{11}$, $B_{11}$, $C_{11}^s$, $D_{11}$, $F_{11}^s$, $H_{11}^s$ and $J_{11}^s$ are stiffness coefficients and are expressed as:

The first variation of work done by axial compressive load $N_0$ can be expressed as:

The governing equations can be derived by introducing Eq. (14) and Eq. (17) into Eq. (13), consequently applying integration by parts, then setting the coefficients $\delta u_0$, $\delta w_0$ and $\delta \varphi$ to zero:

The nonlocal stress state of a specific reference point in a continuum is a function of not only the strain at that exact point but also the strain field throughout the entire body, according to Eringen’s nonlocal elasticity theory. The nonlocal stresses can be derived using Eringen’s differential form \citep{17,18}:

where, $\mu=\left(e_0 a\right)^2$ is the nonlocal parameter, which is responsible for capturing small-scale effects. $\alpha$ and $e_0$ are the internal characteristic length and a material constant, respectively. Literature studies suggest that $e_0 a \leq 2 \mathrm{~nm}$ for single-walled carbon nanotubes \citep{48}. Note that, when $\mu=0$ is considered, the small-scale effect vanishes and the formulation reduces to the classical (local) beam.

The stress resultants are modified by introducing Eq. 10 and (11) into Eq. (19) and subsequent results into Eq. (15):

By replacing Eq. (20) into Eq. (18), the nonlocal governing equations can be expressed in terms of unknown displacements are obtained:

This section presents the Navier solution method to analytically solve the buckling problem of simply supported FGP nanobeams. The simply supported boundary conditions are satisfied as follows:

The displacement components are expressed as trigonometric series, containing unknown coefficients to satisfy the boundary conditions \citep{49}:

where $U_n$, $W_n$, and $\Phi_n$ are the unknown coefficients of displacement. $\beta=n \pi / L$ and $n$ is the mode number.

The nontrivial solution of the buckling of FGP nanobeam is obtained by substituting Eq. (23) into Eq. (21):

The objective of the current proposed hybrid ML approach was to transform the verified solution of the governing formulation of functionally graded porous nanobeam (FGPNB) to a rapid, interpretable surrogate model to assist engineering decisions based on the powerful prediction performance of the hybrid ML approach. In addition, to explore the ability of the hybrid ML approach in predicting complex structural behavior. The main concept of the proposed approach aligns with recent ML implementations in predicting buckling performance of FGP beams and other thin-walled structures \citep{50,51,52}. The novelty of the proposed hybrid ML modeling is combining analytical features, baseline with residual learning, preventing leakage, uncertainty analysis, and explanatory artificial intelligence by SHAP, in addition to multi-objective design within one unified framework. As seen in Figure 2, the data of 540 unique samples were preprocessed, including their $L/h$, PDT, $\alpha_0$, $\mu$, and $\bar{N}$ values. The Python-PyCharm code included standardization and checking the alignment of each sample with its accumulated values, and implemented a one-hot encoder to handle the PDT (UD, FGO, and FGX) as numeric values. These features were the main and basic features used to map the learning approach, then further features were integrated to expose the interaction trends between features and the dominant geometry. These features were generated by creating relations between $L/h$, $\alpha_0$, and $\mu$. However, to minimize the large variation between $L/h$ values, a logarithmic space was used to learn the positive target. However, to avoid using any information from the outer validation fold, all preprocessing procedures were accomplished inside the CV loop. Based on the literature, since ensemble boost models, ANN, and support vector revealed significant prediction performance, the current hybrid ML framework considers implementing all of these models in addition to Gaussian process regression (GPR). The ensemble models included eXtreme Gradient Boosting (XGBoost), histogram gradient boosting (HistGB), random forest (RF), and extra trees. For ANN, the multilayer perceptron (MLP) was used. The support vector regression was implemented based on radial basis (support vector regression-radial basis function, SVR-RBF). However, all algorithms were evaluated two times; first as a direct regression model, then as hybrid ML modeling. Some of the implemented parameters for ML algorithms and their value are illustrated in Table 1. However, the code was prepared and processed based on selecting the best combination of parameters that governed the optimum results.

The hybrid scheme has combined the nonlinear residual-learning model with the regularized Ridge model. Initially, the Ridge regression has estimated the principal trend of the DBL (on a logarithmic scale). Subsequently, the residual was defined as:

where, $t_i$ and $\hat{t}_{\text{base}}$ stand for the logarithmic of the analytical DBL and the prediction of Ridge, respectively. while $x_i$ is equivalent to $x$ deified in Figure 2b.

However, by using the training data only, the parameter of Ridge regularization was selected as:

Following that, the $r_i$, was set as the target to be predicted by the nonlinear model, causing the final hybrid prediction to be redefined as:

where, $\hat{r}_\text{ML}(x)$ defined the predicted residual by ML, while the exponential was used to turn back the converted logarithmic scale to its original scale. Accordingly, Ridge represented the trend of principal response, while the ML corrected the remaining nonlinear discrepancy. However, the best selected residual ML learner model was revealed to be SSVR-RBF, as demonstrated in Section 5. The expression of SVR-RBF was integrated through the LIBSVM library in PyCharm, which used the well-known formulation of SVR network \citep{53}. Nested CV was also used in selecting the model and hyperparameters. The architecture was built by considering five comprised and shuffled folds as an outer layer, which were joined by $L/h$-PDT using 42 random seeds. Each single outer has contained an inner 3-fold search, which created six hyperparameter sampled combination per each model. Consequently, each case has received one-out-of-fold (OOF) prediction from an outer unseen model (untrained model). The OOF was selected because it offers comprehensive evaluation across the entire sample size while preserving validation integrity. Furthermore, the stability of the model was evaluated by repeating 5-fold stratified validation in 10 repetitions, where the leave-one-level-out and learning curve were tested for each $\mu$, $L/h$ and $\alpha_0$, considering three error metrics as R-squared ($R^2$), mean square error (MAE), and root mean square error (RMSE) to measure accuracy. The predicted results were compared with analytical findings, while the uncertainty was implemented using the disjoint split-conformal process, using separate folds for calibration, training, and testing \citep{54}. Furthermore, the SHAP was used to inspect the importance of each feature, while the PDT was combined with SHAP and displayed different categories in different colors. Lastly, as a multi-objective design optimization, the validated hybrid SVR-RBF surrogate model was coupled with NSGA-II. The optimization approach continuously elevated the conservative buckling load $\alpha_0$ to the maximum while considering $\mu = 2$, since the range of $\mu$ was (0, 1, 2, 3, and 4), 2 was selected because it presents the center intermediate value. Besides, to avoid unsupported extrapolation, optimization variables were limited to the analytical domain only. However, because the buckling resistance reduced with increasing porosity, which was concluded by testing several nanobeam designs and recording their porosity, a single design was not able to offer the maximum value of both. Due to this, the Pareto front was used to determine the optimum available candidate to select between larger porosity and greater buckling resistance \citep{55}. The Pareto front was developed based on NSGA-II and the prediction of the trained hybrid model.

In this section, the numerical results are compared with previous results to verify current formulation. Note that, for all studies, $n=1$ was considered.

In Table 2, the DBL values of FGP beams ($\mu=0$) against $\alpha_0$ for different PDTs were compared with the results of Wu et al. \citep{10}. Here, the following material and geometric properties were assumed as: $E=200$ GPa, $v=0.3$ and $h=1$m , slenderness ratio $L/h=10$. For the DBL, the following formulation was used, unless stated otherwise:

In Table 3, DBLs of isotropic nanobeams ($\alpha_0=0$) versus $L/h$ and $\mu$ with the results of Thai and Vo \citep{28} were compared. Furthermore, following formulation was used for DBL:

It can be seen from Table 2 and Table 3, that the obtained results from present theory were in good agreement with previous studies.

The conducted parametric studies to elucidate the effect of nonlocal parameter, slenderness ratio $L/h$, PDT and $\alpha_0$ on DBLs of the FGPNBs is addressed in this section. The following material and geometrical properties were assumed: $E=$200 GPa, $v=$0.3, $L=$10 nm and $b=$1 nm, the first mode number of buckling was considered ($n=$1) and for the DBL, the following formulation was used:

Figure 3 visualizes the variation of DBL values of FGP nanobeam versus $\mu$ for different PDT and $\alpha_0$, and $L/h=$10. It can be seen from Figure 2, DBLs decrease as $\mu$ and $\alpha_0$ values increase for all PDTs.

In Figure 4, the variation of DBLs versus $L/h$ for different $\alpha_0$ is presented for different PDTs and $\mu=2 nm^2$. Besides, it is revealed that an increase in $L/h$ is not only leads to reduction in DBLs but also causes DBL values from distinct $\alpha_0$ converge.

In Figure 5, the variation of DBLs versus $\alpha_0$ and $\mu$ is presented for different PDTs and $L/h=5$. It is seen that, increasing $\alpha_0$ leads to lower DBL values.

In Figure 6, the variation of DBLs versus $L/h$ with and without porosity, for $\alpha_0=0.5$ and $\mu=2 nm^2$ is displayed. Accordingly, as $L/h$ increased the DBLs have decreased as expected. The FGX distribution exhibited the highest DBLs, whereas the FGO type yields the lowest values for FGP nanobeam. This is attributed to the fact that the FGX nanobeam has lower porosity density near the beam surfaces.

This section, present the results obtained by the proposed ML framework based on parametric studies to elucidate the effect of nonlocal parameter, $L/h$, PDT and $\alpha_0$ on DBLs of the FGPNBs is addressed in this section. Figure 7, represent the evaluation of the proposed ML framework. As seen in Figure 7a the best models are listed and evaluated based on RMSE (when RMSE value is closer to zero, it reflects better prediction performance). As seen, spurious performance was recorded by hybrid performance SVR-RBF with RMSE = 0.00729. In overall, hybrid models revealed better performance. To clarify the enhancement which hybrid model deliver, Figure 7b demonstrates the hybrid ablation. It is illustrated that the RMSE value dropped to 0.01298 when direct ML was used, from 0.1017 obtained by the analytical baseline separately. While a further remarkable reduction was accomplished as 0.00729, when hybrid residual learning was used. Indicating enhancement in baseline and direct ML as 93.9\% and 43.9\%, respectively. Figure 7c represent the parity plot for the best model performance (hybrid SVR-RPF). The values of the three implemented error metrices are recorded as $R^2$ = 0.999900, MAE = 0.00309 and RMSE = 0.00729. Like RMSE, when MEA value is closer to absolute zero, better prediction performance is indicated, while $R^2$ approaching 1 indicates better prediction performance. From the trend of scatter points in following the ideal line closely at the entire range of $\bar{N}$, and for all values of $L/h$ it can be proven that the SVR-RBF was not only numerically the best candidate, it also success to capture the complex physics and mechanics across the entire domain. However, the reason of generating parity plot for $L/h$ parameter not any other parameter was due to revealing its significant importance on $\bar{N}$, compared with other parameters based on permutation importance analysis and SHAP analysis. Figure 7d shows the residual plots, as seen in overall at all $L/h$ values the minimum error was recorded by UD PDT, while increasing the complexity of porosity distribution cause to slightly increase the RMSE reflecting slightly lower performance. Yet, in overall the RMSE values for all PDTs are significantly low and accepted under this hybrid ML high validity model using nest CV and OOF. However, at $L/h$ = 5 at all PDTs, indicated to exhibiting the greatest buckling loads and steeper performance gradients.

Figure 8 illustrates the uncertainty and model validation results. As seen in Figure 8a, when the training data increased, the RMSE has reduced from $\approx$ 0.022 to 0.005, indicating superior model performance with high fidelity. Besides, the uncertainty band shows to become thinner, reflecting higher stability. However, the training error preserved its small value through the entire processes. It is seen that there is no variation after training samples greater than $\approx$ 360. By considering the Nested CV for RMSE stability illustrated in Figure 8b, the repeated fixed-model CV and nested CV both recorders very close RMSE values as 0.00501 and 0.00525, respectively. Besides, the error of most folds seems to be $\approx$ 0.01, yet few numbers of higher-error folds are seen. However, it is confirmed from the repeated analysis that the model accuracy does not depend on one favored data partition. In generally efficient stability revealed to be accomplished through the entire range of alternative validation partitions. Figure 8c prove that across the implemented hybrid models there was no noticeable desecrations of the analytical trends for all explored parameters. In contrast, the opposite scenario is seen by various direct models. This illustration confirmed that the hybrid residual learning cause enhancements in mechanical readability and numerical accuracy. The uncertainty interval shown in Figure 8d, shows perfect agreement between prediction and analytical solution through the entire $\bar{N}$ range. A pragmatic coverage of 97.04\% achieved by the 95\% conformal intervals. It is worth to mention that, these intervals only quantified in the current predictive domain and does not validate any untested parameter.

The SHAP analysis was conducted for the best model performance only, therefore Figure 9 demonstrate the SHAP analysis implementations and how each parameter affects the final predicted $\bar{N}$. Figure 9a, shows the importance of each parameter, as seen the most important parameter is $L/h$ while the lowest obtained to be $\mu$. The dominant $L/h$ indicates the robust separation of buckling load across the different four values of $L/h$. The secondary contribution revealed to be accomplished by $\alpha_0$. The importance results by SHAP were exactly same as permutation importance results, signifying the reliability of interpretation. Figure 9b combined SHAP bee swarm analysis with categorizing PDT. As seen a high positive SHAP contribution accomplished by low $L/h$, while higher values decreased the predicted buckling load. On the other hand, a shift to negative region caused by increased $\alpha_0$ and $\mu$. The blue, red and green colors indicate the FGX, FGO and UD, respectively. Figure 9c demonstrate the SHAP interaction maps. The influence of $\alpha_0$ became significantly negative with increasing the value of $\alpha_0$, which also depend on $L/h$. However, larger variation seems to be caused by low $L/h$, indicating stringer coupling or interaction between geometry-porosity at high loads area. As another evidence of the ability of the proposed model, in capturing the coupling effects between parameters rather than one single effect per time, the contribution of $\mu$ with respect to $\alpha_0$ was also explored.

The responses surface illustrated in Figure 10a. The discrete horizontal classes reflect the $L/h$ levels without integrating any extrapolation tool. However, it is seen that increasing $L/h$ and $\alpha_0$, caused a drop in $\bar{N}$. At all PDT classes, the largest resistance to buckling, demonstrate to occur at low $\alpha_0$ and at $L/h$ = 5. In general, structures with FGX shows to have greater predicted resistance compared to FGO, mainly, at large $\alpha_0$. Meanwhile, Figure 10b demonstrates the NSGA-II Pareto front, reflecting the correlation between increasing $\alpha_0$, while preserving good buckling resistance. By considering FGX FGPNB, while maintaining $L/h$ as 5 and $\mu$ as 2, and considering $\alpha_0$ as variable, as proven earlier in analytical solution the resistance to bulking reduced with increasing $\alpha_0$. However, with exploring the optimal scenario, the equal-preference compromise design seems to occur at $\alpha_0$ = 0.499 associated with $\bar{N}_\text{LCB}$ = 1.958, at which a balance between buckling resistance and $\alpha_0$ can be accomplished.

This study developed a hybrid analytical–ML framework for predicting and interpreting the buckling behaviour of FGP nanobeams. The analytical formulation combined exponential shear deformation theory with Eringen’s nonlocal elasticity theory. The governing equations were derived from the minimum total potential energy principle and solved by Navier’s method for simply supported boundary conditions. Three porosity distributions based on the Gibson–Ashby model were considered, and the analytical formulation was verified against published benchmark results. The resulting dataset of 540 configurations was then used to construct a residual-learning surrogate that combined a regularised Ridge baseline with nonlinear ML models. Nested CV, conformal prediction, SHAP analysis, and NSGA-II optimisation were incorporated to evaluate prediction accuracy, quantify uncertainty, interpret parameter effects, and support design exploration.

The principal findings are summarised as follows:

$\bullet$ The nonlocal formulation predicted lower DBLs than the corresponding local limit. The DBL decreased consistently as the nonlocal parameter increased, confirming the importance of small-scale effects in the stability analysis of FGP nanobeams.

$\bullet$ Increasing the porosity coefficient reduced the DBL, although the magnitude of this reduction depended on the porosity distribution. Among the distributions examined, FGX produced the highest buckling resistance, FGO produced the lowest, and UD showed intermediate behaviour.

$\bullet$ The DBL decreased with increasing slenderness ratio. The interaction between the slenderness ratio and porosity coefficient was particularly pronounced in configurations associated with relatively high buckling loads.

$\bullet$ The exponential shear deformation formulation satisfied the traction-free conditions at the upper and lower surfaces without requiring a shear correction factor. It therefore provided a consistent representation of transverse shear deformation within the analytical framework adopted in this study.

$\bullet$ The hybrid SVR-RBF surrogate produced the best predictive performance among the examined models, with an RMSE of 0.00729, an MAE of 0.00309, and an $R^2$ of 0.999900. It also outperformed both the regularised baseline and the corresponding direct ML formulation.

$\bullet$ The learning-curve, repeated-validation, and leave-one-level-out analyses showed that the surrogate retained the principal mechanical trends of the analytical model across the investigated parameter domain. The 95\% conformal prediction intervals achieved an empirical coverage of 97.04\%, indicating reliable uncertainty estimates for configurations within this domain.

$\bullet$ The SHAP analysis identified the slenderness ratio as the most influential input, followed by the porosity coefficient, while the nonlocal parameter had the lowest overall contribution within the selected ranges. These interpretations were consistent with the trends obtained from the analytical parametric analysis.

$\bullet$ Coupling the validated surrogate with NSGA-II produced a Pareto front describing the trade-off between porosity and buckling resistance. For the FGX distribution at $L/h$ = 5 and $\mu$ = 2, the equal-preference solution was obtained at $\alpha_0$ = 0.499, with a predicted DBL of 1.958.

The proposed framework shows how analytical mechanics, residual learning, uncertainty quantification, \sloppy interpretable modelling, and multi-objective optimization can be coordinated within a single design procedure. It provides a computationally efficient means of evaluating FGP nanobeam configurations within the investigated parameter space while retaining the mechanical trends of the underlying analytical formulation.

The present study is limited to the linear static buckling of simply supported FGP nanobeams represented by three prescribed porosity distributions. In addition, the surrogate was trained entirely on analytically generated data and was not independently assessed against experimental measurements or high-fidelity numerical simulations. Future work should therefore consider other boundary conditions, grading patterns, loading conditions, nonlocal formulations, and beam theories. Extending the framework to bending, vibration, thermal buckling, and post-buckling responses, together with validation using finite-element or experimental data, would provide a stronger basis for its application to practical nano structural design.