Selecting the appropriate training technique is a significant step in utilizing intelligent approaches . It becomes even more important when it comes to critical problems like analyzing the bearing capacity of foundations. This study investigates the feasibility of a capable metaheuristic algorithm, called water cycle algorithm (WCA) , for training a multi -layer perceptron (MLP) . The WCA -MLP is applied to a large finite element dataset to predict the settlement. The results of this model are compared with electromagnetic field optimization (EFO) and s huffled complex evolution (SCE) benchmarks. With reference to the obtained Pearson correlation factor s (larger than 0.88 in all stages) , all employed models are suitable for the mentioned objective. Moreover, i t was observed that the training error of the WCA was 5.84 and 3.89% smaller than the EFO and SCE , respectiv ely. Likewise , the accuracy of the WCA -MLP was 1.85 and 2.04% larger in the testing phase. Also, a predictive equation is finally elicited for practical applications in com patible circumstances . Keywords: Bearing capacity ; Settlement measurement ; Artificial neural network; Water cycle algorithm. 1 Introduction With recent advances in computational intelligence, many scholars have replaced traditional methods with economical and accurate machine learning , deep learning [1-7], decision making [8; 9], and artificial intelligence -based tools [10-14]. These novel appro ximation techniques are well employed in various engineering field such as in evaluating the environmental concerns [15-25], implications for natural environmental [26-34], water resources management [35-41], energy efficiency [42-50], structural design [51-61], image processing [62-65], feature selection /extraction [66-70], face recognition [71-74], control performance [75], vibration analysis [76], climate change [77], managing the smart cities [78], project management [79], while in the field of medical science artificial intelligence employed to have a better diagnosis of a particular patients [80-84], early diagnosis of them [85; 86], or medical image classification [87]. There have been many novel algorithm s enhancing the current predictive neural network -based models. Metaheuristic algorithms have been highly regarded in various problems that demand an optimal solution [85]. The hybrid optimization techniques such as differential evolution [88], data -driven robust optimization [89], whale optimization algori thm [90; 91], harris hawks optimization [88; 92], differential edge detection algorithm [93], many -objective sizing optimization [94], fruit fly optimization [95], moth -flame optimization strategy [96; 97], bacterial foraging optimization [98], ant colony optimization [99], particle swarm optimization (PSO) [100-102], chaos enhanced grey wolf opt imization [82], and quantum -enhanced multiobjective large -scale [103]. The efficient determination of ultimate bearing capacity (BC) of foundation s is a significant consideration in designing various structures [104-106]. Khorrami et al. [107] presented a n explicit formulation for ultimate BC for foundations settled on granular soil through an M5’ model tree. As an advantage, the proposed model showed lower uncertainty as well as larger efficiency in comparison with a number of conventional theories. Likew ise, Khorrami and Derakhshani [108] used a hybrid of this model coupled with genetic programming for calculating the BC for a system of footing and cohesionless soil. Sethy et al. [109] introduced adaptive neuro -fuzzy inference system (ANFIS) as a capable appr oximator for computing the BC of rectangular footing with eccentrically loading. Also, in comparison with artificial neural network (ANN), the output pattern of the ANFIS was in a larger agreemen t with the expected one (the coefficients of determination of 0.9024 vs. 0.9118). The competency of support vector machine (SVM) and random forest for estimating the BC of footings placed on rocks was demonstrated by Dutta et al. [110] . Moayedi and Hayati [111] tested several soft computing approaches for explor ing the BC of shallow foundations installed near homogeneous slopes. The findings revealed t he excellence of the feedforward ANN compared to models like tree regression fitting and SVM. In a similar effort, Acharyya and Dey [112] investigated the feasibility of the ANN. They also analyzed the importance of effective factors using this model and found that the angle of internal friction plays the most influential role. Aouadj and Bouafia [113] embedded a new mathematical activation func tion and fed the proposed network by the records of the cone penetration test. Large agreements between the desired and produced outputs, as well as around 30% higher accuracy compared to classical approaches , demonstrated the suitability of the proposed m odel. More applications of the ANNs for this objective can be found in many earlier studies [114-116]. It is well accepted that diverse optimization theories and algorithms can overcome the difficulties that come up with intricate problems [89; 117-127]. Genetic algorithm (GA) is a popular optimization method which was applied by Hamrouni et al. [128] for probabilistic analysis of seismic BC. Saha et al. [129] presented a solution to the ultimate BC problem using symbiosis organisms search. Confirmed by numerical analysis (in the PLAXIS environment), the acceptability of the proposed method was shown for future application s. Jin et al. [130] could accurately study the ultimate BC and critical slip surface of a r ough embedded foundation placed on sands using improved radial movement optimization (IRMO). The efficiency of biogeography -based optimization algorithm (BBO), evolution strategy (ES), and differential algorithm (DE) was investigated by Kashani et al. [131] for the optimal design of the foundation. Gandomi and Kashani [132] professed the superiority of teaching -learning -based optimization algorithm (TLBO) over various swarm -based strategies, such as accelerated PSO, WOA, MFO, for the economical design of shallow foundation s. Moayedi et al. [133] compared two highly popular optimizers of imperialist competition algorithm (ICA) and particle swarm optimization (PSO) for enhancing the performance of the ANN applied to bearing capacity estimation of shallow circular footing. The outputs of the ANNs developed by the PSO and ICA achieved the R2 values of 0.9575 and 0.9467, respectively. Therefore, they concluded the excellence of the PSO -ANN model. A comparison between four metaheuristic strategies of ant colony optimization (ACO), league champion optimization (LCA), whale optimization algorithm (WOA), and moth –flame optimization (MFO) incorporated with ANN in forecasting the stability/failure of a soil -footing system was conducted by Moayedi et al. [134] . Respective accuracies of 94.4, 93.5, 93.9, and 93.9% showed the larger optimization competency of the ACO algorithm. The literature addresses the wide application of metaheuristic algorithm s in the BC calculation [135; 136]. It has been well concluded that these algorithms can efficiently deal with providing optimal solutions, due to their global search capability. Based on the findings of earlier studies concerning the successful use of newly -designed techniques, this study uses water cycle algorithm (WCA) for training an ANN in the BC analysis. The WCA is a robust search strategy that has been broadly employed for different applications [137; 138]. Moreover, the competency of this model is checke d by two benchmark techniques, namely electromagnetic field optimization (EFO) and s huffled complex evolution (SCE) which are known as relatively quicker algorithms . This strategy, i.e., utilizing automatic error mi nimization techniques as the trainer of ANN, has been regarded as a promising solution for complex engineering issues. 2 Methodology 2.1 The WCA technique Most metaheuristic algorithms are nature -inspired, m eaning that the main idea is eli cited from the natural phenomena or the behavior of creatures . The WCA , as is implied by the name, mimicks the way rivers and streams end up with the sea [139]. Melted snow and glaciers flow downhill to form a stream (or a river). Their water is evaporated, becomes clouds, and returns to the earth [140]. The steps required for implementing the WCA can be explained as follows: • Step1: Setting the parameters of the algorithm including 𝐾𝑠𝑟, 𝐾𝑝𝑜𝑝, 𝑑𝑚𝑎𝑥, 𝐼𝑡𝑚𝑎𝑥. • Step1: Scattering the initial population and determining sea, streams, and rivers . Assuming 𝐾𝑝𝑜𝑝 as the total size of the population, 𝐾𝑠𝑟 = 1 + Number of rivers , and 𝐾𝑠𝑡𝑟𝑒𝑎𝑚𝑠 (= 𝐾𝑝𝑜𝑝 - 𝐾𝑠𝑟) as the number of streams, this process is expressed by the below equation : Total population = [ 𝑆𝑒𝑎 𝑅𝑖𝑣𝑒𝑟1 ⋮ 𝑆𝑡𝑟𝑒𝑎𝑚 𝐾𝑠𝑟+1 ⋮ 𝑆𝑡𝑟𝑒𝑎𝑚 𝐾𝑝𝑜𝑝] = [ 𝑥11𝑥21…𝑥𝐾1 𝑥12𝑥22…𝑥𝐾2 ⋮⋮⋮⋮ 𝑥1𝐾𝑝𝑜𝑝𝑥2𝐾𝑝𝑜𝑝…𝑥𝐾𝐾𝑝𝑜𝑝] (1) where river, stream, and the sea represent a particular solution by a 1 × K -dimensional array as [𝑥1, 𝑥2, ..., 𝑥𝐾]. • Step 3: The cost of each member of the existing population is reflected as follows: 𝐶𝑗 = 𝐶𝑜𝑠𝑡𝑗 = f (𝑥1𝑗,𝑥1𝑗,…,𝑥𝐾𝑗) j = 1, 2, …, 𝐾𝑝𝑜𝑝 (2) • Step 4: Given 𝑁𝑆𝑘 as the number of streams discharging in the corresponding rivers or the sea, Equation 3 and 4 give the flow intensity of sea and rivers : 𝐶𝑘 = 𝐶𝑜𝑠𝑡𝑘 - 𝐶𝐹𝐾𝑠𝑟+1 k = 1, 2, …, 𝐾𝑠𝑟 (3) 𝑁𝑆𝑘 = 𝑟𝑜𝑢𝑛𝑑 {|𝐶𝑘 ∑𝐶𝑘𝐾𝑠𝑟 𝑘=1|×𝐾𝑆𝑡𝑟𝑒𝑎𝑚𝑠}, k = 1, 2, …, 𝐾𝑠𝑟 (4) • Step 5: Below relationships describe the streams flow ing into the sea and rivers: 𝑋𝑠𝑡𝑟𝑒𝑎𝑚(𝑡+1) = 𝑋𝑠𝑡𝑟𝑒𝑎𝑚(𝑡) + rand (0, 1) × G × (𝑋𝑠𝑒𝑎(𝑡) - 𝑋𝑠𝑡𝑟𝑒𝑎𝑚(𝑡)) (5) 𝑋𝑠𝑡𝑟𝑒𝑎𝑚(𝑡+1) = 𝑋𝑠𝑡𝑟𝑒𝑎𝑚(𝑡) + rand (0, 1) × G × (𝑋𝑟𝑖𝑣𝑒𝑟(𝑡) - 𝑋𝑠𝑡𝑟𝑒𝑎𝑚(𝑡)) (6) where G is a number varying from 1 and 2 (close to 2). Notably, the streams are allowed to move to the rivers from different directions once C > 1. • Step 6: The movement of the rivers toward downhill (or the sea) can be formulated as follows : 𝑋𝑟𝑖𝑣𝑒𝑟(𝑡+1) = 𝑋𝑟𝑖𝑣𝑒𝑟(𝑡) + rand × G × (𝑋𝑠𝑒𝑎(𝑡) - 𝑋𝑟𝑖𝑣𝑒𝑟(𝑡)) (7) • Step 7: The position of a stream that gives a better -fitted so lution replace s that of the river. • Step 8: Likewise , the position of a river which gives a better -fitted solution replace s that of the sea. • Step 9: Given 𝑑𝑚𝑎𝑥 as a very small value for contro lling the intensification level, BU and B L as the upper and lower bound , respectiv ely, the following procedure checks the conditions of evaporation (for unconstrained problems) : If |𝑋𝑠𝑒𝑎− 𝑋𝑟𝑖𝑣𝑒𝑟𝑗| < 𝑑𝑚𝑎𝑥 or rand < 0.1 j = 1, 2, …, 𝐾𝑠𝑟 – 1 Rain based on Equation 9 End if (8) 𝑋𝑠𝑡𝑟𝑒𝑎𝑚𝑛𝑒𝑤(t + 1) = B L + rand × (BU - BL) (9) As for constrained problems, the WCA uses the below code for enhancing its capability: If |𝑋𝑠𝑒𝑎− 𝑋𝑠𝑡𝑟𝑒𝑎𝑚𝑗| < 𝑑𝑚𝑎𝑥 j = 1, 2, …, 𝑁𝑆𝑘 Rain based on Equation 11 End if (10) 𝑋𝑠𝑡𝑟𝑒𝑎𝑚𝑛𝑒𝑤(t + 1) = X sea + √𝛿 × randn (1, K) (11) where randn is a random number and δ signifies the variance term showing the searching region in the vicinity of the sea . This is worth noting that for preventing premature convergence in such problems, Equation 10 is only implemented for the streams which have direct movements toward the sea. • Step 10: The 𝑑𝑚𝑎𝑥 is decreased as follows: 𝑑𝑚𝑎𝑥(t + 1) = 𝑑𝑚𝑎𝑥(t) - 𝑑𝑚𝑎𝑥(𝑡) 𝐼𝑡𝑚𝑎𝑥 (12) • Step 11: The algorithm is finished if any stopping criteri on is met, otherwise it repeats the process from step 5 [141; 142]. 2.2 Benchmark optimizer Electromagnetic field optimization, as the name indicates, is inspired by the electromagnetic behavior of different poles . In this regard, electromagnets with similar/dissimilar poles repel /attract each other . This algorithm was suggested by Abedinpourshotorban et al. [143] . The i ndividuals of the EFO are electromagnetic particles. The goodness of each particle is the basis of classifying them in to three groups by the name positive field, negative field, and neutral field. The main idea of this technique for improving the solutions is generating and replacing new particles. More clearly, the generated particle with a better fitness magnitude replace s the worst one. For generating a new individual, one particle is randomly chosen from each one of the triple fields . The position and the pole of the neutral particle are first given to the produced individual. The positive and negative electromagnetics also affect it (attraction and repul sion process, respectively ) [144]. More information about the EFO mechanism can be found in the related literature [145; 146]. Shuffled complex evolution was designed by Duan et al. [147] . This algorithm draws on four major concepts including (i) synthesizing probabilistic and deterministic theories, (ii) performing systematic evolution s on so -called containers “ complex es”, (iii) doing competitive evolution s, and (4) shuffling the complex es. Scattering the points is the initial step. Considering the function value of the points, they are then sorted in ascending order. Next, the points are partitioned into a number of complexes where each of them can independe ntly evolve and search the viable space. By applying the modified simplex method of Nelder and Mead [148] for global enhancement , some points are chosen from each unit to form a sub -complex. The larger the fitness of the point is, the more the likelihood of generating offspring is. The worst po ints are replaced by the new offspring [149]. Studies like [150; 151] have explained the SCE in more detail. 3 Data collection As is known, the settlement ( Uy) of a footing is a function of several parameters. In this work, the effect of friction angle (FA), setback distance (SD), elastic modulus (EM), unit weight (UW), dilation angle (DA), applied stress (AS), and Poisson's ratio (PR) is incorporated for measuring the settlement by a series of finit e element analysis executed on a two -layered soil system that bears a shallow footing. Figure 1 depicts the histogram chart of the mentioned parameters. According to the statistical indicators, the values of FA, SD, EM, UW, DA, AS, and PR range in [30, 42] , [1, 7] m, [17500, 65000] 𝑘𝑁 𝑚2, [19.0, 21.1] 𝑘𝑁 𝑚3, [3.4, 11.5], [0.0, 1132.6] 𝑘𝑁 𝑚, and [0.2490, 0.3330]. Also the target parameter (i.e., the Uy) varies from 0 to 0.10 m. (a) (b) (c) (d) (e) (f) (g) (h) Figure 1: Histogram diagram of the parameters of the prepared dataset. With reference to different conditions for this problem, a total of 901 Uys were recorded for the executed stages. The statistical indicators of mean, standard error, standard deviation, and samp le variance are 0.0380, 0.0011, 0.0325, and 0.0011, respectively, for the obtained Uy values. After writing the settlements in front of the corresponding input parameters, the data were permuted randomly and 721 samples (around 80% of records) were selecte d for exploring the Uy behavior. 4 Results and discussion This research is an effort to find an appropriate novel trainer of ANN in examining the bearing capacity of a shallow foundation. To achieve this, a hybrid of WCA and MLP is created and applied to predict the settlement of the foundation. The quality of the prediction is addressed by using three accuracy indicators of Pearson correlation coefficient (R), mean absolute error (MAE), a nd root mean square error (RMSE) . Given 𝑈𝑦 𝑖𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑒𝑑 and 𝑈𝑦 𝑖𝑒𝑥𝑝𝑒𝑐𝑡𝑒𝑑 as the simulated and real Uys, respectively, Equations 13 to 15 formulate the R, MAE, and RMSE. exp expexp 1 exp 122 1( )( ) () -- ()-- predicted ected predicted ectedG y i y y i ypredicted ected i GG y i y y i ypredicted ected iiU U U U R U U U U= ===  (13) exp 11|| ected predictedG y i y i iMAE U U G==− (14) exp2 11[( )] ected predictedG y i y i iRMSE U UG==− (15) 4.1 WCA -MLP results The implementation of any neuro -metaheuristic model consists of several steps. The WCA -MLP is first created by assigning the WCA as the trainer of an MLP network with a 7 × 7 × 1 architecture. As explained, initializing the parameters of the algorithm is the first step. Based on a trial and error process, 𝑑𝑚𝑎𝑥 and 𝐾𝑠𝑟 were set to be 1e -16 and 4, respectively. Also, examining the convergence behavior of the WCA indicated the suitability of 1000 for the maximum number of iterations. The population size (Npop) is another effective factor that was optimized by testing nine values. The converge nce curves of the tested models are shown in Figure 2. In this figure, the objective function (on the y-axis) is the RMSE of the training data. Figure 2: The sensitivity of the WCA performance to the population size . This practice revealed that the best results are y ielded by the WCA -MLP with Npop = 400. The corresponding RMSE was 0.013677. This value, as well as the MAE = 0.0094281 , indicate s an acceptable level of accuracy in understand ing the Uy behavior. As for predicting this pattern, the RMSE and MAE were 0.015175 and 0.010515 , respectiv ely. As is shown in Figure 3, the correlation for both phases is around 90% which implies a good agreement between the expected and modeled Uys. The R indices were exactly 0.90527 and 0.89365 for the training and testing samples. (a) (b) Figure 3: The regression charts of the (a) training and (b) testing results for the WCA -MLP . 4.2 Benchmark results The EFO -MLP and SCE -MLP were developed in the same way and predicted the Uy. The convergence curves of these models are illustrate d in Figure 4. As is seen, the EFO required 5000 iterations to reach a stable convergence during training the MLP. Note that, the Npops for the EFO and SCE were 25 and 10 determined by a trial and error practice. (a) (b) Figure 4: The convergence curve s of the benchmark methods . According to the training re sults, these methods could satisfactorily analyze the relationship between this parameter and related factors. The RMSEs for these two models were close (i.e., 0.014476 and 0.014209 ). Likewise, the MAEs of 0.010383 and 0.010042 indicated a similar level of accuracy for both models in training the ANN. In the testing phase, the RMSEs of 0.015455 and 0.015484 along with the MAEs of 0.011104 and 0.010976 proved that the used models can present a reliable prediction of the intended parameter. Figure 5 depicts the regression charts of these models. The R values of 0.89331 and 0.89743 are obtained for the training data and 0.88973 and 0.89001 are obtained for the testing data of the EFO - MLP and SCE -MLP models . These results demonstrate that the produced Uys are well correlated with the expected values. (a) (b) (c) (d) Figure 5: The regression charts of the training and testing results for the (a and b) EFO -MLP and (c and d) SCE-MLP . 4.3 Comparison From two previous sections, it can be derived that all three metaheuristic algorithms (WCA, EFO, and SCE) could act as a capable trainer for the MLP neural network. In this section, the performance of the proposed algorithm is compared with benchmark ones. Figure 6 shows graphical view s of the training and testing errors calculated for the outputs of the used models. It can be seen that the errors obtained from the WCA -optimized model are more aggregated near the ide al line (Error = 0) . (a) (b) Figure 6: A graphical comparison of the (a) training and (b) testing errors for all used models. -0.07-0.05-0.03-0.010.010.03 0 100 200 300 400 500 600 700Error Data numberEFO-MLP SCE-MLP WCA-MLP -0.07-0.05-0.03-0.010.010.03 0 20 40 60 80 100 120 140 160 180Error Data numberEFO-MLP SCE-MLP WCA-MLP Moreover, in terms of all accuracy ind icators, the prediction of the WCA -MLP was more accurate than EFO -MLP and SCE -MLP. For example, examining the training RMSEs showe d that the error of the EFO -MLP and SCE -MLP is 5.84 and 3.89% larger than WCA -MLP (relative to the WCA - MLP). These values were 1.85 and 2.04 % for the testing data. Conside ring the role of the metaheuristic algorithms in combination with the neural system, it is deduced that the MLP configuration designed by the WCA performs more reliably than the other two search strategies. The time taken by the used algorithms for finding the optimal responses is considered, too. While t he WCA trained the ANN in 4807.8 seconds, the EFO and SCE needed around 48.9 and 546.5 seconds. Notably, the system used for executing the algorithms was a personal computer with the CPU of Intel core i7 with 16 gigs of RAM) . Concerning the reasons for these distinctions, apart from the essen ce of searching strategies , the benchmark models were implemented with simpler configurations (i.e., smaller Npops). 4.4 The formula of the WCA -MLP In this section, t he explicit formula of the WCA -MLP is exhibited as a series of linear/non -linear relationships . Figure 7 shows the architecture of the MLP neural network used for approximating the Uy from FA, SD, EM, UW, DA, AS, and PR . Figure 7: The neural structure of the used predictive model . The Uy is calculated as follows: Uy = -0.023666 × HO 1 + 0.901186 × HO 2 - 0.643029 × HO 3 + 0.945461 × HO 4 - 0.716870 × HO 5 - 0.788679 × HO 6 - 0.104387 × HO 7 + 0.465768 (16) where 𝐻𝑂𝑖= 2 1+ 𝑒−2×𝑠𝑖−1 (17) in which 𝑠1 = 1.084288 × FA + 0.595115 × DA + 0.291941 × UW + 1.247339 × EM + 0.391523 × PR - 0.138593 × SD - 0.272225 × AS - 1.848657 (18) 𝑠2 = -0.249199 × FA - 0.755865 × DA - 0.986040 × UW - 0.972855 × EM - 0.144769 × PR + 0.130731 × SD + 0.909577 × AS + 1.232438 (19) 𝑠3 = 1.023457 × FA - 0.043397 × DA + 0.855424 × UW + 0.984439 × EM + 0.515326 × PR - 0.177398 × SD - 0.608510 × AS - 0.616219 (20) 𝑠4 = -0.158531 × FA 0.271860 × DA - 1.049325 × UW - 1.325108 × EM + 0.375280 × PR - 0.145215 × SD - 0.547336 × AS + 0.000000 (21) 𝑠5 = -0.607473 × FA - 0.701226 × DA - 1.026808 × UW - 0.867166 × EM - 0.420470 × PR - 0.103031 × SD - 0.750380 × AS - 0.616219 (22) 𝑠6 = 0.773196 × FA + 0.180922 × DA - 1.280874 × UW + 0.700894 × EM - 0.298608 × PR + 0.712766 × SD - 0.240574 × AS + 1.23243 8 (23) 𝑠7 = -0.051134 × FA + 0.712179 × DA + 0.311052 × UW - 0.103681 × EM - 0.963517 × PR + 1.000864 × SD - 0.932803 × AS - 1.848657 (24) The numbers used in the above equations represent the weigh ts and biases of the MLP optimally found by the WCA during minimizing the learning error. Equation 16 releases the Uy by doing a linear computation on hidden outputs ( 𝐻𝑂𝑖). As Equation 17 denotes, calculating 𝐻𝑂1, 𝐻𝑂2, ..., and 𝐻𝑂7 consists in obtaining 𝑠1, 𝑠2, ..., and 𝑠7 from Equation s 18 to 24, respectively. In fact, Equation 17 represents a so-called activation function “Tansig” that is used for the neurons in the hidden layer. According to many previous studies, Tansig is a suitable function that can nicely deal with abrupt changes in the dataset [152; 153]. 5 Conclusions The methods based on neural computing have been popularly used for bearing capacity analysis. In this work, the water cycle algorithm, ele ctromagnetic field optimization, and shuffled complex evolution algorithm s were appointed for training an MLP neural network in approximating the settlement of a shallow foundation . The main conclusions are as follows: • Implementing the models by their optimal parameters results in a higher quality of training. • Compared to the EFO and SCE, the WCA needs a considerably larger population for accomplishing the optimization (Npops of 400 vs. 25 and 10 ). Also, the EFO was implemented 5000 times to reach a stable situation. • Based on more than 88% correlation of the outputs in all stages, the used hybrid models can properly capture and reproduce the Uy behavior. • The WCA -MLP surpassed the other two models in terms of all RMSE, MAE, and R accuracy indicators. In other words, the searching strategy of this algorithm could find a more promising solution to the given problem ( i.e., setting the appropriate configuration of the MLP ). • The tested WCA -ANN can be used for accurate analysis of the Uy in practice. Conflict of interest: Authors declare no conflict of interest.