Development of Hybrid Machine Learning Models for Predicting the Critical Buckling Load of I-Shaped Cellular Beams

The principal purpose of this work is to develop three hybrid machine learning (ML) algorithms, namely ANFIS-RCSA, ANFIS-CA, and ANFIS-SFLA which are a combination of adaptive neuro-fuzzy inference system (ANFIS) with metaheuristic optimization techniques such as real-coded simulated annealing (RCSA), cultural algorithm (CA) and shuffled frog leaping algorithm (SFLA), respectively, to predict the critical buckling load of I-shaped cellular steel beams with circular openings. For this purpose, the existing database of buckling tests on I-shaped steel beams were extracted from the available literature and used to generate the datasets for modeling. Eight inputs, considered as independent variables, including the beam length, beam end-opening distance, opening diameter, inter-opening distance, section height, web thickness, flange width, and flange thickness, as well as one output of the critical buckling load of cellular steel beams considered as a dependent variable, were used in the datasets. Three quality assessment criteria, namely correlation coefficient (R), root mean squared error (RMSE) and mean absolute error (MAE) were employed for assessment of three developed hybrid ML models. The obtained results indicate that all three hybrid ML models have a strong ability to predict the buckling load of steel beams with circular openings, but ANFIS-SFLA (R = 0.960, RMSE = 0.040 and MAE = 0.017) exhibits the best effectiveness as compared with other hybrid models. In addition, sensitivity analysis was investigated and compared with linear statistical correlation between inputs and output to validate the importance of input variables in the models. The sensitivity results show that the most influenced variable affecting beam buckling capacity is the beam length, following by the flange width, the flange thickness, and the web thickness, respectively. This study shows that the hybrid ML techniques could help in establishing a robust numerical tool for beam buckling analysis. The proposed methodology is also promising to predict other types of failure, as well as other types of perforated beams.


Introduction
In the field of structural engineering, the use of steel cellular beams has gained widespread attention thanks to their economic and aesthetic advantages [1,2]. Currently, as a result of their particular design, they are widely used in structural systems subject to bending or funicular structures [3][4][5]. The I-shaped cellular steel beams are formed from standard I-section beams by cutting regularly circular openings, a specific geometrical feature that results in numerous advantages, for example, the optimum self-weight-depth ratio, improvement of flexural stiffness, or larger section modulus [3]. Furthermore, in combination with the significant compressive strength of concrete slab, the flexural resistance of composite cellular beam is significantly improved [2,6]. Safety at construction projects is an important matter, which is strongly related to the stability of the constituent structural elements, in other words, the bucking behavior of cellular beams. In the literature, various experimental, analytical, and numerical studies have been performed to investigate the buckling behavior of beams [7][8][9]. As regards to I-shaped beams [5], a numerical investigation was carried out on the buckling modes of cellular beams using ABAQUS software. An important change in the failure loads was observed when subjected to a combination of web distortional and web-post buckling failing modes. In addition, the lateral-torsional buckling behavior of cellular beams was investigated by Sonck and Belis [10] who took into account the effect of the residual stress. In another attempt, the bucking behavior of steel I-shaped cellular beams was discussed and numerous parametric studies were used to assess the variation of the moment-gradient coefficient [4]. This study also provided a modification to evaluate the critical lateral buckling moment of the analyzed beams under different loading conditions. Last but not least, an experimental and analytical study on the behavior and strength of steel beams was performed by Tsavdaridis and D'Mello [11]. An empirical formula with the ability to predict ultimate vertical shear load strength was proposed using a finite element approach in a parametric analysis. In general, the mechanical simulation approach or laboratory experiments can only be applied to a limited case [12]. Moreover, they are costly and require significant amounts of time [13]. Therefore, the development of a straightforward and rapid numerical tool to solve the bucking behavior problem of steel I-shaped cellular beams is imperative [1,12,14].
In the past few decades, machine learning (ML) approaches have been the subject of extensive interest and expanding research in the field of materials science [15], especially in structural engineering [8,9,[16][17][18][19]. Currently, the rapid and constant development of computational software and hardware has facilitated the development of numerous alternative computer-aided data mining approaches. For example, Reddy [20] developed an artificial neural networks (ANN) algorithm for predicting the hardness of concrete reinforced with randomly oriented steel fibers. The proposed ANN model has been compared with selected experiments and has shown an extremely high utilization rate. In addition, other studies [21][22][23][24][25][26] have also confirmed the strong predictability of different ML models for the behavior of structural elements and different materials in the field of mechanics under different solicitations [8,9,[27][28][29][30].
Consequently, the main objective of this paper is the development and comparison of different hybrid ML models to predict the critical buckling load of I-shaped steel cellular beams. In this regard, three hybrid ML models, namely ANFIS-RCSA, ANFIS-CA, and ANFIS-SFLA, which are a combination of adaptive neuro-fuzzy inference system (ANFIS) with metaheuristic optimization techniques such as real-coded simulated Annealing (RCSA), cultural algorithm (CA) and shuffled frog leaping algorithm (SFLA), were developed and applied. Performance of these models was validated and compared using well-known statistical measures such as correlation coefficient (R), root mean squared error (RMSE), and mean absolute error (MAE). Then, the best prediction model was used to perform sensitivity analysis to deduce the factors that greatly affect the predicted critical buckling load.

Description of Cellular Beams and Selection of Variables for Training ML Models
The data used in this work was extracted from a validated finite element (FE) model, previously introduced in the literature by Abambres et al. [31]. The goodness of using data from a validated FE solver instead of experiments for training ML models has been proved in various investigations such as Mallela and Upadhyay [32] and Sadovský and Soares [33]. In this study, a total number of 3645 FE configurations on buckling behavior of I-section beams with circular openings were extracted to form the dataset for training the ML models, involving eight inputs and one output, detailed as follows: The cross section of the beams is I-shaped under the four controlled parameters, section height, flange width, thickness, and web thickness ( Figure 1). As indicated in Table 1, the section height varies from 420 mm to 700 mm (560 mm of mean value and 114.33 mm of standard deviation). The web thickness ranges from 9 mm to 15 mm (mean value of 12 mm and standard deviation of 2.45 mm). The flange width varies from 162 mm to 270 mm (216 mm of mean value and 44.10 mm of standard deviation). The flange thickness ranges from 15 mm to 25 mm (mean value of 20 mm and standard deviation of 4.08 mm).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 22 The data used in this work was extracted from a validated finite element (FE) model, previously introduced in the literature by Abambres et al. [31]. The goodness of using data from a validated FE solver instead of experiments for training ML models has been proved in various investigations such as Mallela and Upadhyay [32] and Sadovský and Soares [33]. In this study, a total number of 3645 FE configurations on buckling behavior of I-section beams with circular openings were extracted to form the dataset for training the ML models, involving eight inputs and one output, detailed as follows: The cross section of the beams is I-shaped under the four controlled parameters, section height, flange width, thickness, and web thickness ( Figure 1). As indicated in Table 1, the section height varies from 420 mm to 700 mm (560 mm of mean value and 114.33 mm of standard deviation). The web thickness ranges from 9 mm to 15 mm (mean value of 12 mm and standard deviation of 2.45 mm). The flange width varies from 162 mm to 270 mm (216 mm of mean value and 44.10 mm of standard deviation). The flange thickness ranges from 15 mm to 25 mm (mean value of 20 mm and standard deviation of 4.08 mm). For a beam with given length, various circular openings had been cut regularly along the beam axis, under three controlled parameters such as beam end-opening distance, opening diameter, and inter-opening distance ( Figure 1). The beam end-opening distance ranges from 12 mm to 718 mm (mean value of 265.36 mm and standard deviation of 157.46 mm). The opening diameter varies from 247 mm to 560 mm (383.56 mm of mean value and 92.98 mm of standard deviation). The interopening distance ranges from 24.70 mm to 274.40 mm (mean value of 112.51 mm and standard deviation of 68.51 mm). Finally, the beam length varies from 4000 mm to 8000 mm (6000 mm of mean value and 1410 mm of standard deviation).  For a beam with given length, various circular openings had been cut regularly along the beam axis, under three controlled parameters such as beam end-opening distance, opening diameter, and inter-opening distance ( Figure 1). The beam end-opening distance ranges from 12 mm to 718 mm (mean value of 265.36 mm and standard deviation of 157.46 mm). The opening diameter varies from 247 mm to 560 mm (383.56 mm of mean value and 92.98 mm of standard deviation). The inter-opening distance ranges from 24.70 mm to 274.40 mm (mean value of 112.51 mm and standard deviation of 68.51 mm). Finally, the beam length varies from 4000 mm to 8000 mm (6000 mm of mean value and 1410 mm of standard deviation).
The cellular beams were subjected to a uniformly distributed load, as schematized in Figure 1. The pinned-pinned boundary condition was applied for the beams. It should be noticed that twisting rotation at both ends of the beams was not allowed to expose the in-plane's rotation and translation influence. The constituent material of the beams is steel under two controlled elastic constants such as the Young's modulus of 210 GPa and the Poisson's coefficient of 0.30. Quad shell elements containing 4 nodes and 6 degrees of freedom at each node (3 for translations and 3 for rotations) were employed to discretize the beams in the FE scheme. The FE simulations and numerical results are obtained using ABAQUS [34]. It has been proved that such element can correctly model thin and moderate thickness shell structures with elastic and elasto-plastic behaviors of materials [35]. The flanges were discretized by 12 elements, whereas the discretization of the web along the span were performed with element size of about H/20 ( Figure 2). In addition, a pinned condition was adopted for discretization of the web-flange connection to form a safe lower bound which was independent of the flange size and type of weld, as recommended by Tsavdaridis and D'Mello [11]. The critical buckling load, or the output of the ML models, varied from 26.40 N/m to 1361.70 N/m (mean value of 225.68 N/m and standard deviation of 80.87 N/m). It should be pointed out that initial imperfections corresponding to geometry or fluctuations of material and residual stresses were not studied [36].   The data used in this work were scaled into the range [0, 1] in order to minimize bias within the dataset between variables. For the jth variable, the expression for scaling is defined as follow: where scaling parameters m and n are also indicated in Table 1. In general, ML related problems often use this technique to transform the dataset into a uniform range to reduce numerical bias [9,37]. It should be noted that a reserve transformation could be deduced for converting data from the scaling space to the raw one. The relationship between input and output variables is plotted in Figure 3. It is observed that no particular correlations are found between two inputs or between inputs and the output.  The joint relationships and histograms between the considered input variables and critical buckling load of the problem (notation is indicated in Table 1).

Adaptive Neuro-Fuzzy Inference System (ANFIS)
AFNIS is a multilayered network based on neural network learning algorithms and fuzzy reasoning [38][39][40][41]. Thanks to its ability to combine the a fuzzy system with numerical power of adaptive neural networks, ANFIS has demonstrated its ability to model many problems in the field of science such as prediction of the geopolymer concrete compressive strength [42], the Marshal parameters of stone matrix asphalt mixtures [43], the buckling damage of structural members [9], the dynamic loadings of power systems [44,45], the buckling damage of steel columns under axial compression [8], the International Roughness Index of pavements [46], the blast-induced air-overpressure in quarry blasting sites [47], the traffic air pollution [48], and the flocculation-dewatering performance of fine mineral tailings [49]. ANFIS has a good ability to learn, build, expand and classify as the advantage of ANFIS is that it allows the fuzzy rules to be extracted from numerical data and is appropriate for formulating a rule base. Moreover, ANFIS can adjust the complex movement of human intelligence into fuzzy systems. However, computation time is one disadvantage of ANFIS, especially in the training process and identification of parameters.

Real-Coded Simulated Annealing (RCSA)
RCSA is a powerful tool for global optimization. On the basis of the similarity between a search algorithm and the process of annealing in metallurgy, the first idea of simulated annealing (SA) appeared in the work of Metropolis et al. [50] as a simulation algorithm. The search algorithm of SA focuses on identifying the solutions without ignoring better solutions that can be found later. Kirkpatrick et al. [51] and Cerny et al. [52] used Metropolis' idea to search for feasible solutions and converge to an optimal solution called "simulated annealing".
Similar to a cooling process, the algorithm simulates a steady temperature decrease until the system converges to a stable state, avoiding the inclusion of defects when cooling too quickly or too slowly [53]. While annealing is the process of first heating of a solid and then cooling it down slowly, in SA, the temperature is kept variable to simulate the heating mechanism. Specifically, the temperature is set high at the initial step and then "cool down" occurs slowly. The initial heating essentially helps to avoid a trap in a local minimum [54]. As the temperature decreases, its new structure becomes increasingly fixed, thus, firmly setting the final properties. In the end, the free energy of the system is minimized, imitating how a minimum is reached during the annealing process, eventually resulting in an optimized solution [55,56]. For more information on SA, the reader is referred to the published studies [57,58]. In this study, RCSA was used to optimize the parameters of the ANFIS to develop a hybrid model, namely ANFIS-RCSA.

Cultural Algorithm (CA)
Derived from observations of cultural evolution in real life, CA is a computational and optimization algorithm [59,60]. It is helpful to deal with complex calculations or optimization of the searching process. CA is inspired from the human culture development at two levels of evolutionary. The first one, the macro level, takes place in the beliefs space and the second one, the micro level, occurs in the space of people [61].
In brief, CA consists of three components. First, it is the evolved social population and the evaluation, reproduction, and modification mechanisms. Secondly, the belief space represented for the bias was acquired by the population during the problem-solving process. The communication protocol, the third component, is where the interaction of the population and their beliefs is identified [62].
Precisely, the process starts by initializing, firstly, the belief, as well as the population spaces. In the next step, the algorithm is repeated for each generation before achieving the termination condition. The performance function is used to evaluate the individuals. The influence function and the acceptance function are mediums to communicate between two levels of CA. The two functions identify and select the individuals from the current population to influence the belief space. The experiences of selected ones are generalized and applied to adjust the current beliefs via the updated function. The generated beliefs are further used to guide and influence the next generation through the evolutionary process [60].
CA is popularly used to perform complex calculations of the optimization search algorithm [63] and in many other applications [62,64]. In this study, CA was used to optimize the parameters of the ANFIS to develop a hybrid model, namely ANFIS-CA.

Shuffled Frog Leaping Algorithm (SFLA)
SFLA belongs to the metaheuristic algorithms based on frog's behavior. It was first introduced in the work of Eusuff and Lansey [65] to optimize the water distribution network. The frog's behavior in seeking food is the main concept of SFLA. Moreover, the advantages of other optimization methods, such as particle swarm optimization (PSO) or genetic algorithm (GA), could also be found in SFLA algorithm [66,67]. In the SFLA optimization process, the global and local search mechanisms are properly compromised. Several subsets (memeplex) divide the frog population and perform in the local search. In each memeplex, the worst frog is replaced by the best one [68,69].
Numerous complex optimization problems can be solved by SFLA, such as multimodal problems, nonlinear and non-differential [67,[70][71][72]. This method can achieve better results through the distribution of water resources, the fuzzy or multivariable PID controller, and the repair of bridge-deck [73,74]. In this study, SFLA was used to optimize the parameters of the ANFIS to develop a hybrid model, namely ANFIS-SFLA.

Performance Indicators
The performance of three hybrid ML models was evaluated by three statistical criteria, namely RMSE, R, and MAE [7,[75][76][77]. The R-value indicates the correlation between the predicted values, given by the ML models, and the actual values [78]. RMSE is the average squared difference between the predicted and actual values [19,79,80]. MAE measures the average of absolute difference between outputs (predicted values) and targets (actual values) [81][82][83]. On the one hand, higher R indicates a better model [84]. On the other hand, lower RMSE and MAE indicate a better model [85,86]. The formulas are defined as follow [87][88][89][90]: where n refers to the number of samples in the dataset; p i and q are the values and means of the predicted critical buckling load, respectively; and v i and v are the individual and mean values of the actual critical buckling load, respectively.

Methodology Flow Chart
The work's methodology ( Figure 4) is described below, containing five main steps: Dataset preparation to train ML models: The input and output parameters were used to create a complete set of data. A number of 70% data (2551 training samples) were extracted from the initial dataset for training the ML models. The remaining 30% data (1094 testing samples) were used for validation the AI models.

3.
Training models: The ML models were trained using the training dataset. Three ML algorithms based on ANFIS with three optimization methods RCSA, CA, and SFLA. The concepts of these models were introduced in the previous sections. This step was repeated until the models successfully trained within a preselected tolerance error criterion.

4.
Models validation: After successfully training three ML models, the validation process was performed using the testing dataset. The models were verified using different statistical measures such as RMSE, MAE. and R.

5.
Sensitivity analysis: After validation of the ML models, the sensitivity of input parameters was performed using the best model to identify the influence between input parameters and the critical buckling load of I-shaped cellular beams.

Building the Hybrid ML Models
Performance of a ML model highly depends on the choice of its parameters. In this study, the parameters for the ANFIS model was selected and used to develop the hybrid models after trial and error test on training dataset (Table 2). In addition, the performance of the model can be affected by the sampling strategy of datasets; thus, the sampling strategy was analyzed by randomly picking data points in the whole dataset to reduce the effect of randomness on the model's performance. Over 100 times of simulation was ran to get the best performance of the models corresponding to the highest values of R and lowest values of MAE and RMSE. Different optimization techniques such as RCSA, CA, and SFLA were applied for optimizing the parameters of ANFIS and the corresponding parameters of these techniques are summarized in Tables 3-5, respectively. Evaluation of cost functions with respect to R and MAE using ANFIS-RCSA, ANFIS-CA, and ANFIS-SFLA for predicting the critical buckling load of I-shaped beams are presented in Figure 5. For each model, the stopping condition was chosen as 1000 iterations. It is shown that ANFIS-RCSA reaches the state of

Building the Hybrid ML Models
Performance of a ML model highly depends on the choice of its parameters. In this study, the parameters for the ANFIS model was selected and used to develop the hybrid models after trial and error test on training dataset (Table 2). In addition, the performance of the model can be affected by the sampling strategy of datasets; thus, the sampling strategy was analyzed by randomly picking data points in the whole dataset to reduce the effect of randomness on the model's performance. Over 100 times of simulation was ran to get the best performance of the models corresponding to the highest values of R and lowest values of MAE and RMSE. Different optimization techniques such as RCSA, CA, and SFLA were applied for optimizing the parameters of ANFIS and the corresponding parameters of these techniques are summarized in Tables 3-5, respectively. Evaluation of cost functions with respect to R and MAE using ANFIS-RCSA, ANFIS-CA, and ANFIS-SFLA for predicting the critical buckling load of I-shaped beams are presented in Figure 5. For each model, the stopping condition was chosen as 1000 iterations. It is shown that ANFIS-RCSA reaches the state of convergence over 100 iterations, whereas it is about 800 iterations for ANFIS-CA and ANFIS-SFLA.

Validating the Hybrid ML Models
Validation of the hybrid models was done using both training and testing datasets, as shown in Figure 6. The ideal regression lines are in discontinuous, whereas the R-values between predicted and actual critical buckling load are given in Table 6. The results are obtained after 1000 iterations for all hybrid ML models. The clouds of data points are given in color in order to better illustrate the concentration of points. The correlation results analysis shows that the R values vary from 0.927 to 0.962 for the training dataset and from 0.920 to 0.960 for the testing dataset. These results display a very good prediction performance of all the models, but ANFIS-SFLA is found to be the best predictor compared with other models.
A detailed summary of the prediction capability of ANFIS-RCSA, ANFIS-CA, and ANFIS-SFLA is presented in Table 6, whereas "Error Mean" is the mean of all relative errors, "Error StD" is the standard deviation of relative errors and "Slope" represents the slope between the ideal regression line and the linear fit given by the models. The results also confirm that ANFIS-SFLA is superior to ANFIS-RCSA and ANFIS-CA in predicting the critical buckling load of I-shaped beams.

Lowest concentration
Highest concentration On the basis of the training phase analysis, it is observed that selection of optimal parameters for training the predictive ML models is highly important before any further investigation. The stopping iteration of a model is also an important factor as it could directly affect the computation time, as well as the efficiency of an algorithm. A more appropriate convergence criterion should be included in the ANFIS-SFLA code in order to better track the information of the cost function (i.e., R, RMSE, and MAE). By default, the 1000 iterations criterion as the stopping condition is selected because it provided a reasonable computation time. Such criterion should be modified in further studies such that the error between the cost function of iteration (i + 1) and iteration (i) exposes no more than a given threshold (i.e., 1% for example). In general, it can be concluded that ANFIS-SFLA is the best ML predictor, following by ANFIS-CA and ANFIS-RCSA in the training phase.

Validating the Hybrid ML Models
Validation of the hybrid models was done using both training and testing datasets, as shown in Figure 6. The ideal regression lines are in discontinuous, whereas the R-values between predicted and actual critical buckling load are given in Table 6. The results are obtained after 1000 iterations for all hybrid ML models. The clouds of data points are given in color in order to better illustrate the concentration of points. The correlation results analysis shows that the R values vary from 0.927 to 0.962 for the training dataset and from 0.920 to 0.960 for the testing dataset. These results display a very good prediction performance of all the models, but ANFIS-SFLA is found to be the best predictor compared with other models. and actual critical buckling load are given in Table 6. The results are obtained after 1000 iterations for all hybrid ML models. The clouds of data points are given in color in order to better illustrate the concentration of points. The correlation results analysis shows that the R values vary from 0.927 to 0.962 for the training dataset and from 0.920 to 0.960 for the testing dataset. These results display a very good prediction performance of all the models, but ANFIS-SFLA is found to be the best predictor compared with other models.
A detailed summary of the prediction capability of ANFIS-RCSA, ANFIS-CA, and ANFIS-SFLA is presented in Table 6, whereas "Error Mean" is the mean of all relative errors, "Error StD" is the standard deviation of relative errors and "Slope" represents the slope between the ideal regression line and the linear fit given by the models. The results also confirm that ANFIS-SFLA is superior to ANFIS-RCSA and ANFIS-CA in predicting the critical buckling load of I-shaped beams.

Highest concentration
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 22    A detailed summary of the prediction capability of ANFIS-RCSA, ANFIS-CA, and ANFIS-SFLA is presented in Table 6, whereas "Error Mean" is the mean of all relative errors, "Error StD" is the standard deviation of relative errors and "Slope" represents the slope between the ideal regression line and the linear fit given by the models. The results also confirm that ANFIS-SFLA is superior to ANFIS-RCSA and ANFIS-CA in predicting the critical buckling load of I-shaped beams.
Histograms of the scaled critical buckling load of I-shaped beams of both actual and predicted values issued from ANFIS-RCSA, ANFIS-CA, and ANFIS-SFLA algorithms is presented in Figure 7 for both training and testing datasets The results show that a high concentration of the values is found in the [0, 0.5] interval, which emphasize the values around 0.1. This is also confirmed by the red color from Figure 6. It is worth noticed that the better correlation between actual and predicted critical buckling load is represented by the height of bars. An equally height of bars for a given value of critical buckling load means a perfect correlation (R = 1). More precisely, the error histograms between actual and predicted critical buckling load simulated by ANFIS-RCSA, ANFIS-CA, and ANFIS-SFLA are presented in Figure 8. The errors are found in the [−0.1, 0.1] range, which speaks for about 10% of error. ANFIS-SFLA is again proved as the best predictor, following by ANFIS-CA and ANFIS-RCSA, respectively. This could be easily observed by the highest value centered on 0 (ANFIS-SFLA) and the narrow curves (ANFIS-SFLA and ANFIS-CA) of about 2% of error. In general, in terms of six error criteria, three models perform well, but ANFIS-SFLA is shown to be superior to ANFIS-RCSA and ANFIS-CA for predicting the critical buckling load of I-shaped beams. On the basis of a comparison with the previous published study by Abambres et al. [31] using the same datasets but different ML model, namely ANN, it can be observed that although the R values of ANFIS-SFLA model in this study (0.96) is smaller than those of ANN model (0.9999), the results of ANFIS-SFLA model are more reasonable for predicting the critical buckling load of I-shaped beams as the R value of 0.9999 of ANN model in the previous study indicated the serious overfitting problem which should be avoided in predictive problems. Moreover, the advantage of the ANFIS-SFLA model is that it used an optimization technique, namely SFLA, to optimize the parameters of ANFIS model for avoiding the overfitting and improve the accuracy of the predictive model [65]. In addition, with respect to the computational time for modeling, the training phase of the ANFIS-SFLA model took about 45 min using 64 bit Windows computer with processor Intel(R) Core (TM) i5-3210M CPU @ 2.50GHz, whereas only a fraction of second was consumed for evaluating new beam configuration. Once constructed, the ANFIS-optimized model can simulate new configurations without generating new beam configuration, as needed in classical FE modeling.
the overfitting and improve the accuracy of the predictive model [65]. In addition, with respect to the computational time for modeling, the training phase of the ANFIS-SFLA model took about 45 min using 64 bit Windows computer with processor Intel(R) Core (TM) i5-3210M CPU @ 2.50GHz, whereas only a fraction of second was consumed for evaluating new beam configuration. Once constructed, the ANFIS-optimized model can simulate new configurations without generating new beam configuration, as needed in classical FE modeling.   an optimization technique, namely SFLA, to optimize the parameters of ANFIS model for avoiding the overfitting and improve the accuracy of the predictive model [65]. In addition, with respect to the computational time for modeling, the training phase of the ANFIS-SFLA model took about 45 min using 64 bit Windows computer with processor Intel(R) Core (TM) i5-3210M CPU @ 2.50GHz, whereas only a fraction of second was consumed for evaluating new beam configuration. Once constructed, the ANFIS-optimized model can simulate new configurations without generating new beam configuration, as needed in classical FE modeling.

Sensitivity Analysis
In this section, the sensitivity analysis was performed in favor of quantifying the influence of these inputs to the predicted critical buckling load using the best-found model, namely ANFIS-SFLA. The procedure of the sensitivity analysis is detailed in the following: A new input space of eight dimensions was constructed based on the probability density function of each variable for evaluating the sensitivity analysis. The new input space composes of only eleven values for each input variable ranging from its smallest to highest quantile with a resolution of 10 (for instance, [0, 10,20,30,40,50,60,70,80,90,100]). When calculating the sensitivity for a considered input variable, the values were taken successively from the eleven previously mentioned quantiles, whereas other inputs remained at their median. Therefore, the method obtains information on the deviation of output solution in varying statistically the considered input variable. In this work, the reference case refers to the one where all input variables stayed at their 50% quantile level. For illustration purpose, the more changes on the output solution as compared with the reference case, the more influence the considered input is. Deviation in output solution when evaluating the sensitivity of the jth input variable can be written as the following [48]: where Y re f is the output of the reference case, Y j i is the output solution using jth input at its ith quantile. Figure 9 shows the evaluation of deviation for eight input variables used in this study. Obviously, all considered input variables govern the buckling capacity of the cellular beam. The critical buckling Appl. Sci. 2019, 9,5458 14 of 21 capacity of the beam decreases when the length, L, and the opening diameter, D, increases, but it increases when the inter-opening distance d, section height H, web thickness t web , flange width w flange , and flange thickness t flange increase. In addition, the beam end-opening distance, d 0 , does not significantly affect the buckling capacity of the beam. However, it could be seen that the deviation in output solution is nonlinear in function of input variables, especially L, t web , t flange , and w flange . In case of the length L, a high order of amplitude is observed as compared with other input variables.
where ref Y is the output of the reference case, j i Y is the output solution using jth input at its ith quantile. Figure 9 shows the evaluation of deviation for eight input variables used in this study. Obviously, all considered input variables govern the buckling capacity of the cellular beam. The critical buckling capacity of the beam decreases when the length, L, and the opening diameter, D, increases, but it increases when the inter-opening distance d, section height H, web thickness tweb, flange width wflange, and flange thickness tflange increase. In addition, the beam end-opening distance, d0, does not significantly affect the buckling capacity of the beam. However, it could be seen that the deviation in output solution is nonlinear in function of input variables, especially L, tweb, tflange, and wflange. In case of the length L, a high order of amplitude is observed as compared with other input variables. Using the deviation previously presented, the final influence level of each input is calculated by summing all deviations based on the following expression: Figure 9. Evaluation of deviation in output θ (%) for eight input variables, respectively.
Using the deviation previously presented, the final influence level of each input is calculated by summing all deviations based on the following expression: The influence level of all input variables obtained by the ANFIS-SFLA model was scaled into the range of [0, 1] as plotted in Figure 10b. It is seen that the most influenced variable is the beam length L (Θ L = 0.429), following by the flange width w flange (Θ w flange = 0.215), the flange thickness t flange (Θ t flange = 0.130), the web thickness t web (Θ t web = 0.121), the inter-opening distance d (Θ d = 0.047), the opening diameter D (Θ D = 0.046), the section height H (Θ H = 0.008), and the beam end-opening distance d 0 (Θ d 0 = 0.003), respectively.
In parallel, linear statistical correlation between each input and output was also calculated and compared with the results of sensitivity analysis obtained by using ANFIS-SFLA (Figure 10a). It is found that the most dominant relationship with the output is the beam length variable (R = 0.667), followed by the flange width (R = 0.375), web thickness (R = 0.332), flange thickness (R = 0.209), opening diameter (R = 0.107), section height (R = 0.092), inter-opening distance (R = 0.072), and beam end-opening distance (R = 0.023), respectively.
It can be observed that there is a good agreement between the results of the ANFIS-SFLA algorithm with the classical regression approach ( Table 7). The classification order of importance obtained by ANFIS-SFLA is also indicated in Table 7. On the basis of the results from ANFIS-SFLA and linear statistical correlation, it can be concluded that the beam length highly affects the critical buckling load of circular openings I-shaped beam, followed by the flange width, flange thickness, and web thickness, respectively. and linear statistical correlation, it can be concluded that the beam length highly affects the critical buckling load of circular openings I-shaped beam, followed by the flange width, flange thickness, and web thickness, respectively.
Overall, without solving complex equations, an approach using hybrid ML models could also indicate these highly nonlinear relationships in an efficient manner and avoid too much computational cost.   Overall, without solving complex equations, an approach using hybrid ML models could also indicate these highly nonlinear relationships in an efficient manner and avoid too much computational cost.

Conclusions
Accurately predicting the critical buckling load is crucial in the field of structural engineering applications. In this study, different hyrbid ML models, namely ANFIS-RCSA, ANFIS-CA, and ANFIS-SFLA, i.e., a combination of ANFIS with metaheuristic optimization techniques such as RCSA, CA and SFLA, respectively, were used to predict the critical buckling load of I-shaped cellular steel beams with circular openings. A total of 3645 FE configurations on buckling behavior of I-section beams with circular openings were extracted to generate the training and testing datasets. Out of these, a total of eight input variables, namely the beam length, flange width, flange thickness, web thickness, inter-opening distance, opening diameter, section height, and beam end-opening distance, as well as one output variable of the critical buckling load of I-shaped cellular steel beams were used in the modeling. Various validation criteria, namely R, RMSE, and MAE were used to validate and compare the models. In addition, a sensitivity analysis was also done to evaluate the importance of input variables for predicting the critical buckling load of I-shaped cellular steel beams using the ML models.
The results showed that the three hybrid proposed models performed well but ANFIS-SFLA is the best predictor as compared with other models (ANFIS-CA and ANFIS-SFLA). The sensitivity results showed that the length of beam is a dominant factor for predicting the critical buckling load of I-shaped cellular steel beams using the ML models, followed by the flange width, flange thickness, web thickness, inter-opening distance, opening diameter, section height, and the beam end-opening distance, respectively.
Overall, the hybrid ML technique developed in this study offers a powerful tool to accurately predict the critical buckling load of I-shaped cellular steel beams. However, the application of the ML model in engineering practice and design is not often an appropriate way. In further research, the current ANFIS-SFLA prediction function could be exploited to derive an explicit empirical equation for use in practical engineering design and analysis. For this purpose, the sensitivity analysis previously presented could play an important role in deriving the desired formula by searching, first, the empirical coefficients between the output and the most relevant inputs.

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