# Prediction of the Coefficient of Pressure Fluctuations during the Hydraulic Jump Using ELM, GMDH, and M5MT

Authors: Tzu-Chia Chen; Biju Theruvil Sayed; Maria Jade Catalan Opulencia; Raed H. C. Alfilh; Maki Mahdi Abdulhasan; Sayed Hashmat Sadat
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

Pressure fluctuations are a critical phenomenon that can endanger the safety and stability of hydraulic structures, especially stilling basins. Hence, the accurate estimation of the dimensionless coefficient of pressure fluctuations (CP′) is critical for hydraulic engineers. This study proposed predictive soft computing models to estimate CP′ on sloping channels. Therefore, three robust soft computing methods, including extreme learning machine (ELM), group method data of handling (GMDH), and M5 model tree (M5MT), were used to estimate CP′. The results revealed that ELM was more accurate than GMDH and M5MT methods when comparing statistical indices, including correlation coefficient (CC), root mean square error (RMSE), mean absolute error (MAE), scatter index (SI), index agreement (Ia), and BIAS values. The performance of ELM was found to be more accurate (CC = 0.9183, RMSE = 0.0067, MAE = 0.0051, SI = 11.88%, Ia = 0.9569) when compared with the results of GMDH (CC = 0.8818, RMSE = 0.0078, MAE = 0.0058, SI = 13.89%, Ia = 0.9361) and M5MT (CC = 0.6883, RMSE = 0.0120, MAE = 0.0090, SI = 21.28%, Ia = 0.7905) in the testing stage. In addition, the BIAS values revealed that ELM slightly overestimated the values of CP′, especially at the peak point compared with GMDH and M5MT results. Overall, the suggested soft computing techniques worked well for predicting pressure fluctuation changes in the hydraulic jump.


---

## Body

## 1. Introduction

Stilling basins are the most widely used dissipation hydraulic structures of large dams. Energy dissipation in stilling basins is related to hydraulic jumps with high turbulent flow. Hydraulic jump is the common phenomenon for flow energy dissipation in the stilling basins. This phenomenon transforms the supercritical flow into a subcritical flow at a short distance. This is also accompanied by large-scale turbulence, surface waves, air entrainment into the flow, an increase in flow depth, and considerable energy dissipation in the water flow [1, 2]. The turbulent flow within stilling basins is related to the movement of large-scale vortices and severe pressure fluctuations, which may cause significant damage in stilling basins through cavitation, erosion, lifting force, and material fatigue [3]. The importance of pressure fluctuations in hydraulic jumps was revealed after the destruction of the Karnafoli and Malpaso dams in Bangladesh and Mexico, respectively [4]. Therefore, pressure fluctuations in hydraulic structures have included a considerable volume of hydraulic engineers’ investigations. Hydraulic models are standard tools for measurements of pressure fluctuations. Experimental studies were carried out for the analysis of characteristics of pressure fluctuations at the flip buckets and stilling basins.Soft computing methods are alternative tools compared with traditional regression approaches and analytical solutions for modeling complicated problems and systems [5]. Moreover, soft computing techniques are much more cost-effective and have less run time than experimental studies. It is noteworthy that soft computing models are newer tools than experimental studies. The capability of different soft computing approaches has been proven for modeling CP′ in the hydraulic jump [6–8]. These models were developed and assessed using experimental data. The suggested soft computing methods were more accurate than conventional regression schemes for predicting CP′ during the hydraulic jump. Moreover, Samadi et al. [9] modeled dynamic pressure distribution in flip bucket spillways using a variety of soft computing approaches and provided simple, high-precision mathematical expressions for estimating it.According to the literature review’s assessment, there were few investigations into the performance of soft computing techniques for estimating pressure fluctuations in hydraulic structures. In addition, the applications of ELM, GMDH, and M5MT methods have not yet been reported in estimating pressure fluctuations. However, the results mentioned above show that soft computing methods are a desirable alternative to model pressure fluctuations in hydraulic structures.The ELM algorithm is an advanced method that is constructed on a single hidden layer feed-forward neural network (SFLN) [10]. The ELM improves training time and accuracy by transforming training data into fixed-length batches and only updating the weight without retraining the trained samples. The ELM approach has been successfully reported for longitudinal dispersion coefficients in water pipelines [11], compressive strength concrete estimation [12], and predicting total dissolved gas [13].The group method of data handling is a set of induction techniques that can be used to make mathematical modeling of multi-parametric datasets. GMDH algorithm uses an inductive approach to sort and order more complex polynomial models, with an external criterion selecting the best result. The application of GMDH has been widely and successfully carried out in scour depth estimation and water-related engineering problems [14–16].M5MT is a common decision tree method that uses piecewise multiple linear regression equations to approximate nonlinear functions. The M5MT was applied to predict drought events [17], scour depth [15, 18], alga growth in reservoirs [19], the capacity of shallow foundations [20], and groundwater modeling [21].The major goal of this study is to investigate the usefulness and capability of ELM, GMDH, and M5MT for predictingCP′ is taking place during the hydraulic jump on sloping channels. The author’s knowledge indicated that ELM, GMDH, and M5MT applications had not been investigated yet for prediction CP′. Therefore, the current research in the first study evaluates the ability of ELM, GMDH, and M5MT to predict CP′ occurring in hydraulic jumps on sloping channels.


## 2. Materials and Methods

The data description and soft computing methods, including ELM, GMDH, and M5MT, are briefly explained in the following subsections.


### 2.1. Data Description

Data were collected from reliable and published experimental results of Gunal [22]. During the hydraulic jump phenomenon in the laboratory, pressure fluctuations occurred on the sloping channel. Gunal’s [22] experiments were conducted in a 91 cm wide and 320 cm long rectangular flume. The inclined angles of the sloping channel were set to 10, 20, and 30 degrees. The cross-sectional representation of the hydraulic jump downstream of the sloping channel is displayed in Figure 1.Figure 1
The cross-sectional representation of a hydraulic jump downstream of a sloping channel.Guven [7] determined the relationship between P¯′2 as root mean square value of the pressure fluctuations with characteristics of the flow and geometry of channels as follows:(1)P¯′2=fyt,y1,V1,x,g,ρw,θ,where P′2¯ is calculated from the following equation:(2)P′2¯=1N∑i=1NPx,y,n∆t−P¯x,y2,where P¯x,y is the mean pressure, Px,y,n∆t is the value instantaneous of pressure, N is the number of data collected in a discrete-time series, and ∆t is the sampling time step.In (1), yt is the depth of tail-water downstream of the stilling basin; y1 is the gate opening of sluice gate; V1 is the velocity issuing from the gate; x is the distance between the horizontal and inclined parts’ intersection; g is the acceleration of gravity; ρw is water’s mass density; and θ is the angle of the sloping section of the flume. Finally, the dimensional analysis provided the dimensionless functional form of (1) as follows [7]:(3)CP′=fFr1,xy1,yty1,θ,where CP′=P′2¯/V12/2g=RMS/V12/2g and Fr1=V1/gy1 are the dimensionless coefficient of pressure fluctuations and inflow Froude number, respectively. CP′ is an important factor used to describe pressure fluctuations.


### 2.2. Data Preparation for the Development of the Soft Computing Approaches

Overall, 112 values were obtained forCP′ under different flow conditions by recording pressure data. The statistical indices in the testing and training data sets should be similar; otherwise, we will be unable to test our models under certain hydrodynamic conditions that are designated for the train set [23]. Therefore 80% of all data is randomly chosen for training, and the rest is used to assess the created models. The major statistical parameters of the testing and training data set are shown in Table 1. As seen in Table 1, the statistical parameters for testing and training subsets were close.Table 1
The main statistical parameters of training and testing data set for developing the proposed soft computing methods.
DatasetVariablesParameterMinMaxAvgSt. dev.TrainingInputFr14.9448.6626.7121.473x/y10.75031.87511.9937.761yt/y16.75011.1118.9921.721θ102020.3338.670OutputCP′0.0270.0960.0560.017TestingInputFr14.9448.6626.3251.357x/y10.75028.33312.5508.731yt/y16.75011.1118.5091.613θ103018.6368.888OutputCP′0.0280.0900.0570.018The variables ofFr1, x/y1, yt/y1, and θ are considered as input variables of soft computing methods for the prediction of CP′.


### 2.3. Extreme Learning Machine (ELM) and Model Development

Although the standard backpropagation neural network has several advantages, it also has some limitations, including tuning the parameter settings by selecting the number of hidden layers, momentum coefficient, and learning rate values. It suffers from the iterative learning process used to determine the weights, which takes an extended period. In addition, there is a possibility of a local minimum. By contrast, ELM provides it with the advantages of rapid convergence, fewer parameters to tweak, and a high degree of generalization [24].The ELM comprises a single hidden layer that containsL nodes. The following equation is valid for N arbitrarily separate samples xi,ti, with an activation function gx, randomly distributed weights (wi), randomly distributed biases (bi), and output weights (βi).(4)∑i=1Lβigwixj+bi=Oj,j=1,….N,where β=β1,……,βIT is a vector containing the output weights of the hidden layer of L nodes and the output node. It is possible to approximate N samples with zero error using typical SLFN with L hidden nodes and activation function gx. This means that(5)∑j=1NOj−tj=0.In other words, there areβi, wi, and bi such that(6)∑i=1Lβigwixj+bi=tj,j=1,….N.Let(7)H=gwixi+b1⋯gwLxi+bL⋮⋱⋮gwixN+b1⋯gwLxN+bL,β=β1T…βLTT.And(8)T=t1T…tLTT.Then,(9)Hβ=T.H is referred to as the neural network’s hidden layer output matrix. Huang et al. [24] prove L ≤ N hidden neurons as necessary for infinitely differentiable activation functions. H remains unchanged once the biases and weights of hidden nodes are established, and only β needs to be estimated. The least-squares solution with a minimum norm for (9) is as follows:(10)minHw1wL,b1bLβ−T.H is square if the neurons in the hidden layer equal the training set (N), and β may be calculated by inverting H. Nevertheless, to obtain higher generalization, the number of hidden nodes is modified, and it may be less than N. Then, the Moore–Penrose generalized inverse of H must be utilized in this case.(11)β^=H†T,where H† can be used to validate the algorithm because it is the Moore–Penrose generalized inverse of H. The ELM model that has been constructed has the following simple general form for the prediction of CP′:(12)CP′=11+expInW×InV+BHIT×OutW,where InV is the input variables, InW is the input weight matrix, BHI is the hidden neuron vector’s bias, and OutW is the output weight vector. The following relations were obtained for the ELM model to calculate CP′.(13)InV=Fr1yty1xy1θ,InW=−0.41−0.170.88−0.61−0.56−0.98−0.920.020.96−0.24−0.65−0.94−0.81−0.91−0.66−0.140.440.680.550.160.61−0.88−0.93−0.86−0.31−0.65−0.020.360.920.980.490.740.57−0.320.53−0.90−0.070.06−0.39−0.86,BHI=−0.680.29−0.990.550.30−0.840.72−0.24−0.97−0.96,OutW=−0.050.03−0.01−0.090.02−0.02−0.06−0.010.040.01.


### 2.4. Group Method of Data Handling (GMDH) and Model Development

The GMDH polynomial neural network introduced by Ivakhnenko [25] is a feed-forward neural network. This algorithm is a self-organizing system that is used to search progressively for the optimal solution to complicated nonlinear problems. This algorithm approximates the relationship between input and output based on quadratic polynomials. Therefore, GMDH creates new neurons in each layer by connecting pairs of neurons with quadratic polynomials. The goal of a mathematical model problem of GMDH is to discover a function (f^) that can be used to approximate an original function (f) in order to estimate the model’s output (y^) for an assumed input vector X including n input variables [26]. For this, given n data instances of multi-input single-output data pairs, the following results are obtained:(14)yi=fX)=f(xi1,xi2,xi3,…,xin;i=1,2,…,M.A mathematical formulation describes the general equation between input-output variables. The goal now is to design a GMDH network in such a way that the square of the deviation output and the estimated output is as little as possible, which means:(15)∑i=1M’f^xi1,xi2,xi3,…,xin−yi2⟶min.The mathematical expression between input and output variables can be described by a sophisticated discrete variant of the Volterra function called the Kolmogorov–Gabor polynomial. This series is presented in the following forms.(16)yw0+∑i=1nwixi+∑i=1n∑j=1nwijxixj+∑i=1n∑j=1n∑k=1nwijkxixjxk+⋯.It is employed in this study to calculate the GMDH network’s quadratic polynomial, which can be represented as [27]:(17)y^=Gxi,xj=w0+w1xi+w2xj+w3xixj+w4xi2+w5xj2.The mathematical expression (18) shows how neurons in a GMDH network are linked together to create the equation between input-output variables. Using the least-squares regression method, the weighting coefficient values of (17) are determined. This means that the deviation output, y, and the one that is calculated, y^, for each pair of xi and xj input variables are as minimal as possible. Thus, the weighting coefficients values of the quadratic function Gi are determined in order to optimize the fit of the output throughout the entire number of sample data pairs, that is [27].(18)E=∑i=1M’yi−Gi2M⟶min.In its ordinary form, the GMDH method considers all possible combinations of two input variables from a total ofn input variables in order to generate the regression polynomial in the form of (17) that fits the dependent data yi,i=1,2,…,M the best in a least-squares scheme. As a consequence, the initial layer of the GMDH network’s architecture will be selected with Cn2=nn−1/2 input neurons for the creation of the quadratic polynomial based on observations yi,xip,xiq;i=1,2,…,M for varied p,q∈1,2,…,n. To put it another way, it is currently possible to generate M data triples yi,xip,xiq;i=1,2,…,M from observations by utilizing such p,q∈1,2,…,n in the following way:(19)x1px1qy1x2px2qy2xmpxmqym.By replacing the quadratic sub-expression in the shape of (17) for each row of M data triples, the matrix equation shown below is simple to construct.(20)AW=Y,where W and Y are the vectors of unknown quadratic polynomial weighting coefficients in (17) and a vector containing the observed values of outputs, respectively.(21)W=w0,w1,w2,w3,w4,w5T,Y=y1,y2,y3,…,yMT.The superscriptT denotes the matrix’s transposition. It is self-evident that:(22)A=1x1px1qx1px1qx1p2x1q21x2px2qx2px2qx2p2x2q21xmpxmqxmpxmqxmp2xmq2.To solve regression analysis problems, we employ the least-squares approach.(23)W=ATA−1ATY.This returns the vector with the best quadratic (17) weighting coefficients for the entire set of M data triples. It is worth noting that this technique is replicated for each neuron in the subsequent hidden layer, depending on the network’s connectivity structure. Using the GMDH algorithm, the following GMDH network was obtained for the prediction CP′ (see Figure 2).Figure 2
The GMDH network for predictionCP′.With respect to Figure2, the following equations were obtained for the prediction CP′ using the GMDH approach.The first layer output:(24)CP′∗11=0.1101−0.0110yty1+3.6904θ+6.2904×10−4yty12+8.3932×10−5θ2−3.0321×10−4yty1θ,CP′∗31=0.1253−0.0180Fr1−5.5504×10−4θ+0.0013Fr12+9.1284×10−5θ2−3.02542×10−4Fr1θ,CP′∗41=0.0979−0.0079Fr1+0.0014xy1+8.1446×10−5Fr12−6.6223×10−5xy12+5.4575×10−5Fr1xy1,CP′∗61=0.0109+0.01502yty1+0.0011xy1−0.0011yty12−6.2790×10−5xy12+5.8831×10−5yty1xy1.The second layer output:(25)CP′∗12=0.1382−2.6925CP′∗41−1.7543CP′∗61−1.3848×103CP′∗412−1.3189×103CP′∗612+2.7568×103CP′∗41CP′∗61,CP′∗220.1671−0.7202CP′∗11−4.8455CP′∗61+31.2869CP′∗112+70.9947CP′∗612−41.3749CP′∗11CP′∗61,CP′∗320.1845−1.1091CP′∗31−5.1238CP′∗61+36.6985CP′∗312+75.5852CP′∗612−45.3275CP′∗31CP′∗61.The output layer:(26)CP′∗130.0337−0.4779CP′∗12+0.1849CP′∗22+6.7083CP′∗122+5.1880CP′∗222−0.2282CP′∗12CP′∗22.Finally, the value ofCP′ is calculated from the above formulations of GMDH. As seen with replacing the input variables with the GMDH formulation, the CP′∗13 approximate the outcome GMDH network for prediction of CP′.


### 2.5. M5 Model Tree (M5MT) and Model Development

M5MT is one of the most widely used decision tree techniques in data-driven modeling. In the M5MT, the entire input domain is recursively partitioned into subdomains, with each subdomain being predicted using a multiple linear regression model. The graphical M5MT is constructed from a root node, the number of binary branches, a group of inner nodes (splits), and a number of terminal nodes (leaves) [28]. For this reason, the resulting tree model has a clear decision structure and is understandable for everyone.The constructing, pruning, and smoothing of the tree are the three main components of the M5MT algorithm. The splitting criterion is used to construct the primary tree. The expected reduction in error resulting from evaluating each attribute at the node is calculated using this splitting criterion. A measure of the error at a node is defined as the standard deviation of the class values that reach the node. After that, the attribute with the highest anticipated error reduction is chosen. The following formula is used to get the standard deviation reduction (SDR) for M5MT:(27)SDR=sdT−∑iTiT×sdTi,where T denotes the set of instances that reach the node, Ti is the result of separating the node according to the attribute chosen, and sd denotes the standard deviation [29]. To stop splitting, either only a few examples remain, or their class values are less than 5% of the initial instance set’s standard deviation. It is possible to encounter an over-fitting problem based on training data during the creation tree process. Pruning is a technique that has been employed in trying to alleviate this difficulty in the past. It decreases the size of the model tree by deleting splits that don’t meaningfully enhance prediction.For the test data, the pruning algorithm uses an estimate of the predicted error at each node. To begin, the absolute deviation between the observed and estimated output values for each of the training cases entering the node is averaged. Due to the fact that the trees were constructed specifically for this dataset, the average will underrate the predicted error for new instances. This is compensated for by multiplying the output value by the factor (n + v)/(n − v). where n denotes the size of the training instances received at the node [30]. In addition, v is the number of model attributes that signify the output value at that node. The leaf node can be omitted if the estimated error is lower at the parent. As a result, this multiplication is done to make sure that new data, rather than training data, don’t get underestimated by the multiplication.Quinlan [31] describes a smoothing approach that employs the leaf model to determine the estimated value during the smoothing phase. The value is then smoothed by merging it with the linear model’s estimated value for each node along the path back to the root. This requires the following calculation:(28)P′=np+kqn+k,P′ is a prediction that exceeds the upper node; p is a prediction from below that is passed to the current node; q is the model’s predicted value at the node; n represents the number of training instances that have made it to the preceding node, and Wang and Witten [28] constant is denoted by k. Using the M5MT created, the graphical tree, is shown in Figure 3.Figure 3
M5MT created the regression tree to predictCP′.The decision rules concerned with Figure3 are as follows:(29)LM4=CP′=−0.0041×Fr1+0.0792.As seen, four linear rules were obtained for the predictionCP′. Concerning the values of Fr1, the appropriate rule was selected then CP′ was calculated.


## 2.1. Data Description

Data were collected from reliable and published experimental results of Gunal [22]. During the hydraulic jump phenomenon in the laboratory, pressure fluctuations occurred on the sloping channel. Gunal’s [22] experiments were conducted in a 91 cm wide and 320 cm long rectangular flume. The inclined angles of the sloping channel were set to 10, 20, and 30 degrees. The cross-sectional representation of the hydraulic jump downstream of the sloping channel is displayed in Figure 1.Figure 1
The cross-sectional representation of a hydraulic jump downstream of a sloping channel.Guven [7] determined the relationship between P¯′2 as root mean square value of the pressure fluctuations with characteristics of the flow and geometry of channels as follows:(1)P¯′2=fyt,y1,V1,x,g,ρw,θ,where P′2¯ is calculated from the following equation:(2)P′2¯=1N∑i=1NPx,y,n∆t−P¯x,y2,where P¯x,y is the mean pressure, Px,y,n∆t is the value instantaneous of pressure, N is the number of data collected in a discrete-time series, and ∆t is the sampling time step.In (1), yt is the depth of tail-water downstream of the stilling basin; y1 is the gate opening of sluice gate; V1 is the velocity issuing from the gate; x is the distance between the horizontal and inclined parts’ intersection; g is the acceleration of gravity; ρw is water’s mass density; and θ is the angle of the sloping section of the flume. Finally, the dimensional analysis provided the dimensionless functional form of (1) as follows [7]:(3)CP′=fFr1,xy1,yty1,θ,where CP′=P′2¯/V12/2g=RMS/V12/2g and Fr1=V1/gy1 are the dimensionless coefficient of pressure fluctuations and inflow Froude number, respectively. CP′ is an important factor used to describe pressure fluctuations.


## 2.2. Data Preparation for the Development of the Soft Computing Approaches

Overall, 112 values were obtained forCP′ under different flow conditions by recording pressure data. The statistical indices in the testing and training data sets should be similar; otherwise, we will be unable to test our models under certain hydrodynamic conditions that are designated for the train set [23]. Therefore 80% of all data is randomly chosen for training, and the rest is used to assess the created models. The major statistical parameters of the testing and training data set are shown in Table 1. As seen in Table 1, the statistical parameters for testing and training subsets were close.Table 1
The main statistical parameters of training and testing data set for developing the proposed soft computing methods.
DatasetVariablesParameterMinMaxAvgSt. dev.TrainingInputFr14.9448.6626.7121.473x/y10.75031.87511.9937.761yt/y16.75011.1118.9921.721θ102020.3338.670OutputCP′0.0270.0960.0560.017TestingInputFr14.9448.6626.3251.357x/y10.75028.33312.5508.731yt/y16.75011.1118.5091.613θ103018.6368.888OutputCP′0.0280.0900.0570.018The variables ofFr1, x/y1, yt/y1, and θ are considered as input variables of soft computing methods for the prediction of CP′.


## 2.3. Extreme Learning Machine (ELM) and Model Development

Although the standard backpropagation neural network has several advantages, it also has some limitations, including tuning the parameter settings by selecting the number of hidden layers, momentum coefficient, and learning rate values. It suffers from the iterative learning process used to determine the weights, which takes an extended period. In addition, there is a possibility of a local minimum. By contrast, ELM provides it with the advantages of rapid convergence, fewer parameters to tweak, and a high degree of generalization [24].The ELM comprises a single hidden layer that containsL nodes. The following equation is valid for N arbitrarily separate samples xi,ti, with an activation function gx, randomly distributed weights (wi), randomly distributed biases (bi), and output weights (βi).(4)∑i=1Lβigwixj+bi=Oj,j=1,….N,where β=β1,……,βIT is a vector containing the output weights of the hidden layer of L nodes and the output node. It is possible to approximate N samples with zero error using typical SLFN with L hidden nodes and activation function gx. This means that(5)∑j=1NOj−tj=0.In other words, there areβi, wi, and bi such that(6)∑i=1Lβigwixj+bi=tj,j=1,….N.Let(7)H=gwixi+b1⋯gwLxi+bL⋮⋱⋮gwixN+b1⋯gwLxN+bL,β=β1T…βLTT.And(8)T=t1T…tLTT.Then,(9)Hβ=T.H is referred to as the neural network’s hidden layer output matrix. Huang et al. [24] prove L ≤ N hidden neurons as necessary for infinitely differentiable activation functions. H remains unchanged once the biases and weights of hidden nodes are established, and only β needs to be estimated. The least-squares solution with a minimum norm for (9) is as follows:(10)minHw1wL,b1bLβ−T.H is square if the neurons in the hidden layer equal the training set (N), and β may be calculated by inverting H. Nevertheless, to obtain higher generalization, the number of hidden nodes is modified, and it may be less than N. Then, the Moore–Penrose generalized inverse of H must be utilized in this case.(11)β^=H†T,where H† can be used to validate the algorithm because it is the Moore–Penrose generalized inverse of H. The ELM model that has been constructed has the following simple general form for the prediction of CP′:(12)CP′=11+expInW×InV+BHIT×OutW,where InV is the input variables, InW is the input weight matrix, BHI is the hidden neuron vector’s bias, and OutW is the output weight vector. The following relations were obtained for the ELM model to calculate CP′.(13)InV=Fr1yty1xy1θ,InW=−0.41−0.170.88−0.61−0.56−0.98−0.920.020.96−0.24−0.65−0.94−0.81−0.91−0.66−0.140.440.680.550.160.61−0.88−0.93−0.86−0.31−0.65−0.020.360.920.980.490.740.57−0.320.53−0.90−0.070.06−0.39−0.86,BHI=−0.680.29−0.990.550.30−0.840.72−0.24−0.97−0.96,OutW=−0.050.03−0.01−0.090.02−0.02−0.06−0.010.040.01.


## 2.4. Group Method of Data Handling (GMDH) and Model Development

The GMDH polynomial neural network introduced by Ivakhnenko [25] is a feed-forward neural network. This algorithm is a self-organizing system that is used to search progressively for the optimal solution to complicated nonlinear problems. This algorithm approximates the relationship between input and output based on quadratic polynomials. Therefore, GMDH creates new neurons in each layer by connecting pairs of neurons with quadratic polynomials. The goal of a mathematical model problem of GMDH is to discover a function (f^) that can be used to approximate an original function (f) in order to estimate the model’s output (y^) for an assumed input vector X including n input variables [26]. For this, given n data instances of multi-input single-output data pairs, the following results are obtained:(14)yi=fX)=f(xi1,xi2,xi3,…,xin;i=1,2,…,M.A mathematical formulation describes the general equation between input-output variables. The goal now is to design a GMDH network in such a way that the square of the deviation output and the estimated output is as little as possible, which means:(15)∑i=1M’f^xi1,xi2,xi3,…,xin−yi2⟶min.The mathematical expression between input and output variables can be described by a sophisticated discrete variant of the Volterra function called the Kolmogorov–Gabor polynomial. This series is presented in the following forms.(16)yw0+∑i=1nwixi+∑i=1n∑j=1nwijxixj+∑i=1n∑j=1n∑k=1nwijkxixjxk+⋯.It is employed in this study to calculate the GMDH network’s quadratic polynomial, which can be represented as [27]:(17)y^=Gxi,xj=w0+w1xi+w2xj+w3xixj+w4xi2+w5xj2.The mathematical expression (18) shows how neurons in a GMDH network are linked together to create the equation between input-output variables. Using the least-squares regression method, the weighting coefficient values of (17) are determined. This means that the deviation output, y, and the one that is calculated, y^, for each pair of xi and xj input variables are as minimal as possible. Thus, the weighting coefficients values of the quadratic function Gi are determined in order to optimize the fit of the output throughout the entire number of sample data pairs, that is [27].(18)E=∑i=1M’yi−Gi2M⟶min.In its ordinary form, the GMDH method considers all possible combinations of two input variables from a total ofn input variables in order to generate the regression polynomial in the form of (17) that fits the dependent data yi,i=1,2,…,M the best in a least-squares scheme. As a consequence, the initial layer of the GMDH network’s architecture will be selected with Cn2=nn−1/2 input neurons for the creation of the quadratic polynomial based on observations yi,xip,xiq;i=1,2,…,M for varied p,q∈1,2,…,n. To put it another way, it is currently possible to generate M data triples yi,xip,xiq;i=1,2,…,M from observations by utilizing such p,q∈1,2,…,n in the following way:(19)x1px1qy1x2px2qy2xmpxmqym.By replacing the quadratic sub-expression in the shape of (17) for each row of M data triples, the matrix equation shown below is simple to construct.(20)AW=Y,where W and Y are the vectors of unknown quadratic polynomial weighting coefficients in (17) and a vector containing the observed values of outputs, respectively.(21)W=w0,w1,w2,w3,w4,w5T,Y=y1,y2,y3,…,yMT.The superscriptT denotes the matrix’s transposition. It is self-evident that:(22)A=1x1px1qx1px1qx1p2x1q21x2px2qx2px2qx2p2x2q21xmpxmqxmpxmqxmp2xmq2.To solve regression analysis problems, we employ the least-squares approach.(23)W=ATA−1ATY.This returns the vector with the best quadratic (17) weighting coefficients for the entire set of M data triples. It is worth noting that this technique is replicated for each neuron in the subsequent hidden layer, depending on the network’s connectivity structure. Using the GMDH algorithm, the following GMDH network was obtained for the prediction CP′ (see Figure 2).Figure 2
The GMDH network for predictionCP′.With respect to Figure2, the following equations were obtained for the prediction CP′ using the GMDH approach.The first layer output:(24)CP′∗11=0.1101−0.0110yty1+3.6904θ+6.2904×10−4yty12+8.3932×10−5θ2−3.0321×10−4yty1θ,CP′∗31=0.1253−0.0180Fr1−5.5504×10−4θ+0.0013Fr12+9.1284×10−5θ2−3.02542×10−4Fr1θ,CP′∗41=0.0979−0.0079Fr1+0.0014xy1+8.1446×10−5Fr12−6.6223×10−5xy12+5.4575×10−5Fr1xy1,CP′∗61=0.0109+0.01502yty1+0.0011xy1−0.0011yty12−6.2790×10−5xy12+5.8831×10−5yty1xy1.The second layer output:(25)CP′∗12=0.1382−2.6925CP′∗41−1.7543CP′∗61−1.3848×103CP′∗412−1.3189×103CP′∗612+2.7568×103CP′∗41CP′∗61,CP′∗220.1671−0.7202CP′∗11−4.8455CP′∗61+31.2869CP′∗112+70.9947CP′∗612−41.3749CP′∗11CP′∗61,CP′∗320.1845−1.1091CP′∗31−5.1238CP′∗61+36.6985CP′∗312+75.5852CP′∗612−45.3275CP′∗31CP′∗61.The output layer:(26)CP′∗130.0337−0.4779CP′∗12+0.1849CP′∗22+6.7083CP′∗122+5.1880CP′∗222−0.2282CP′∗12CP′∗22.Finally, the value ofCP′ is calculated from the above formulations of GMDH. As seen with replacing the input variables with the GMDH formulation, the CP′∗13 approximate the outcome GMDH network for prediction of CP′.


## 2.5. M5 Model Tree (M5MT) and Model Development

M5MT is one of the most widely used decision tree techniques in data-driven modeling. In the M5MT, the entire input domain is recursively partitioned into subdomains, with each subdomain being predicted using a multiple linear regression model. The graphical M5MT is constructed from a root node, the number of binary branches, a group of inner nodes (splits), and a number of terminal nodes (leaves) [28]. For this reason, the resulting tree model has a clear decision structure and is understandable for everyone.The constructing, pruning, and smoothing of the tree are the three main components of the M5MT algorithm. The splitting criterion is used to construct the primary tree. The expected reduction in error resulting from evaluating each attribute at the node is calculated using this splitting criterion. A measure of the error at a node is defined as the standard deviation of the class values that reach the node. After that, the attribute with the highest anticipated error reduction is chosen. The following formula is used to get the standard deviation reduction (SDR) for M5MT:(27)SDR=sdT−∑iTiT×sdTi,where T denotes the set of instances that reach the node, Ti is the result of separating the node according to the attribute chosen, and sd denotes the standard deviation [29]. To stop splitting, either only a few examples remain, or their class values are less than 5% of the initial instance set’s standard deviation. It is possible to encounter an over-fitting problem based on training data during the creation tree process. Pruning is a technique that has been employed in trying to alleviate this difficulty in the past. It decreases the size of the model tree by deleting splits that don’t meaningfully enhance prediction.For the test data, the pruning algorithm uses an estimate of the predicted error at each node. To begin, the absolute deviation between the observed and estimated output values for each of the training cases entering the node is averaged. Due to the fact that the trees were constructed specifically for this dataset, the average will underrate the predicted error for new instances. This is compensated for by multiplying the output value by the factor (n + v)/(n − v). where n denotes the size of the training instances received at the node [30]. In addition, v is the number of model attributes that signify the output value at that node. The leaf node can be omitted if the estimated error is lower at the parent. As a result, this multiplication is done to make sure that new data, rather than training data, don’t get underestimated by the multiplication.Quinlan [31] describes a smoothing approach that employs the leaf model to determine the estimated value during the smoothing phase. The value is then smoothed by merging it with the linear model’s estimated value for each node along the path back to the root. This requires the following calculation:(28)P′=np+kqn+k,P′ is a prediction that exceeds the upper node; p is a prediction from below that is passed to the current node; q is the model’s predicted value at the node; n represents the number of training instances that have made it to the preceding node, and Wang and Witten [28] constant is denoted by k. Using the M5MT created, the graphical tree, is shown in Figure 3.Figure 3
M5MT created the regression tree to predictCP′.The decision rules concerned with Figure3 are as follows:(29)LM4=CP′=−0.0041×Fr1+0.0792.As seen, four linear rules were obtained for the predictionCP′. Concerning the values of Fr1, the appropriate rule was selected then CP′ was calculated.


## 3. Results and Discussion

Statistical indices such as correlation coefficient (CC), root mean square error (RMSE), mean absolute error (MAE), index of agreementIa, scatter index (SI), and BIAS index was employed to assess the qualitative evaluation of the developed suggested models.(30)CC=∑i=1nxi−x¯yi−y¯∑i=1n’xi−x¯2∑i=1n’yi−y¯2,RMSE=∑i=1N’xi−yi2N,MAE=1N∑i=1Nxi−yi,Ia=1−∑i=1n’xi−yi2∑i=1n’xi−x¯+yi−y¯2,SI=RMSEx¯×100,BIAS=y¯−x¯,where xi and yi denote measured and predicted values, x¯ and y¯ represent the average of measured and predicted values, and N denotes the number of the dataset. In Table 2, the statistical index values for the proposed soft computing methods in the training and testing stages for the prediction of CP′.Table 2
The statistical indices for ELM, GMDH, and M5MT for the prediction ofCP′.
ApproachData setCCRMSEMAEIaSI (%)BIASELMTraining0.97700.00460.00400.98418.000.0008ELMTesting0.91830.00670.00510.956911.880.0020GMDHTraining0.91350.00740.00660.951212.85−0.0005GMDHTesting0.88180.00780.00580.936113.890.0016M5MTTraining0.83090.01080.00900.846718.81−0.0003M5MTTesting0.68830.01200.00900.790521.280.0014Table2 shows that the ELM model performed CP′ prediction with the lowest errors (RMSE, MAE, and SI) and higher coefficient correlation (CC) and index agreement (Ia) than the GMDH and M5MT models in both the training and testing stages. As observed in Table 2, the ELM method with CC = 0.9183, RMSE = 0.0067, and MAE = 0.0051 has the best prediction accuracy compared with GMDH and M5MT. On the other hand, it can be deduced from Table 2 that GMDH (CC = 0.8818, RMSE = 0.0078, and MAE = 0.0058) and M5MT (CC = 0.6883, RMSE = 0.0120, and MAE = 0.0090) have the second and third level of accuracy for the prediction CP′ in the testing stage. Furthermore, the scatter index (SI = 11.88%) for the ELM method is smaller than the SI for GMDH (SI = 13.89%) and M5MT (SI = 21.28%). It is noteworthy that the BIAS values of ELM for training and testing data are 0.0008 and 0.0020, respectively, indicating a slightly overestimated CP′. However, the BIAS values of GMDH in training and testing were −0.0005 and 0.0016, and for M5MT, they were −0.0003 and 0.0014 for the training and testing stages, respectively. Therefore, it can be concluded that, overall, ELM was more conservative for prediction CP′ than GMDH and M5MT.By comparing two soft computing approaches, including the ELM (as the best model) and the M5MT (as the worst model) in the testing stage, it can be shown that the ELM produced significantly lower errors in predictionCP′ than the M5MT. The values of RMSE and MAE obtained via ELM showed that they decreased by about 79.10% and 76.47% compared with the M5MT in the testing stage, respectively.Besides, the CC and Ia values for ELM increased by 33.42% and 21.05% compared with M5MT. The error values of RMSE and MAE from GMDH in the testing stage showed that these values improved by about 53.85% and 55.17% compared with those obtained from M5MT, respectively. Moreover, the CC and Ia values of GMDH increased by about 28.11% and 18.42% compared with M5MT. Based on error measures, it is clear that the ELM model was better than the GMDH and M5MT in terms of accuracy.The mathematical expressions for the predictionCP′ were formulated from the proposed soft computing models. It seems the M5MT result was more transparent and easier to use than ELM and GMDH. M5MT provided four multivariate linear models for the prediction of CP′. Although the derived rules from M5 had less accuracy than ELM and GMDH, these rules are more straightforward for predicting CP′. On the other hand, GMDH presented complex mathematical formulations for the prediction of CP′. The GMDH used two hidden layers of neurons for the prediction of CP′. The GMDH utilized the quadratic polynomial of input variables and the combination of the best neurons in each layer to predict CP′. ELM generated the coefficients matrixes for the prediction of CP′. It seems the ELM mathematical shapes were simple compared with the GMDH equations.M5MT method provided four simple equations for the prediction ofCP′ based on dividing the domain of the problem. It is worth mentioning that M5MT selected the appropriate rule concerning only one input variable (i.e., Froude number) and used two variables, including Fr1 and x/y1 in their rules. However, the ELM and GMDH used all the independent variables to generate the predictive expressions for estimating CP′ values. Although M5MT has led to the creation of simple rules for predicting CP′, these simple formulas cause M5MT to be less accurate than ELM and GMDH expressions. M5MT divided the domain of input variables into four subdomains and represented four regression equations to estimate CP′. In fact, separating the domain of the problem into the local subdomains and combining their results caused an improvement in the accuracy compared with the single equation.Figures4–9 compare observed values of CP′ versus estimated values obtained from ELM, GMDH, and M5MT during the training and testing stages. As observed in scatter plots, the predicted CP′ values from ELM were more concentrated around the ideal line (the 45-degree line).Figure 4
Scatter plots of measured and estimated values ofCP′ for ELM during (a) training stage and (b) testing stage.
(a)(b)Figure 5
Variation ofCP′ with ELM during the training stage (a) and testing stage (b).
(a)(b)Figure 6
Scatter plots of measured and estimated values ofCP′ for GMDH during (a) training stage and (b) testing stage.
(a)(b)Figure 7
Variation ofCP′ with GMDH during the training stage (a) and testing stage (b).
(a)(b)Figure 8
Scatter plots of measured and estimated values ofCP′ for M5MT during (a) training stage and (b) testing stage.
(a)(b)Figure 9
Variation ofCP′ with M5MT during the training stage (a) and testing stage (b).
(a)(b)In addition, the variation ofCP′ values obtained by soft computing methods versus observed CP′ demonstrated the capability of estimating CP′ by the proposed approaches. These figures graphically confirmed the higher ELM accuracy than GMDH and M5MT.These figures indicated that ELM performed better in the training and testing stages than GMDH and M5MT. As shown, the ELM reasonably estimatedCP′ in the training and testing stages (bias = 0.0008 and 0.0020). The remarkable point is that in Figures 4 and 5, the peak values of CP′ by ELM were estimated well. In contrast, GMDH and M5MT slightly underestimated the values of CP′, especially at the peak CP′ values. From the comparison between Figures 4–9, it can be deduced that ELM is more skillful and accurate than GMDH and M5MT in the prediction of CP′.The present study results were compared with earlier research conducted by Samadi et al. [8]. They used three data-driven algorithms, such as MARS, GEP, and CART, for the prediction of CP′. Samadi et al. [8] indicated that CART results for the prediction of CP′ have three and four nonterminal and terminal nodes, respectively. The CART tree structure used only the Fr1 variable with threshold values of 6.31, 5.18, and 5.65 for the CART tree structure. The present study used another common decision tree technique, namely M5MT, to predict CP′.It is worth noting that M5MT and CART are two common decision tree algorithms used for regression problems. In the present study, M5MT, among the four input variables included,Fr1, x/y1, yt/y1, and θ, only selected the Fr1 variable for prediction of CP′. As illustrated in Figure 3, M5MT’s tree structure uses three nonterminal and four terminal nodes. In addition, the splitting values for M5MT were 6.31, 5.18, and 5.65. A comparison of the tree structures of the decision tree (i.e., M5MT) and the proposed CART presented by Samadi et al. [8] revealed that the two models have similar structures. They employed a splitting variable (Fr1) and the same splitting values for the generation of regression trees. As a result, the graphical structures of these two decision tree algorithms for selecting the independent variable (Fr1) and creating four if-then rules are similar. However, M5MT provided linear regression functions while CART presented constant values in their terminal nodes. This was the main distinctive characteristic of the differences between M5MT and CART concerning the nature of CART and M5MT. It should be noted that due to the nature of the M5MT algorithm, the Fr1 parameter was an attribute that caused error reduction for the prediction of CP′ in the regression tree obtained by M5MT. In addition, the present study’s findings about selecting the Fr1 parameter is completely consistent with the results of the CART tree provided by Samadi et al. [8].Also, the statistical measures indicated the results of M5MT and CART were more or less the same. M5MT with RMSE = 0.0120 and MAE = 0.0090, compared with CART [12] with RMSE = 0.012, MAE = 0.009, was similar results. However, M5 generated multivariate linear regression equations while CART generated the constant values for the prediction of CP′. The linear equations by M5MT caused the flexibility and generalizability of M5 to improve the prediction of CP′ compared with the constant values of CP′ that was yielded by CART. Compared with GMDH in the present study with GEP [12], both algorithms provided nonlinear mathematical expressions for the prediction of CP′. The accuracy of GMDH with RMSE = 0.0078 and MAE = 0.0058 is slightly better than GEP results with RMSE = 0.008 and MAE = 0.006. Moreover, compared with the formulas provided by GEP and GMDH, GEP requires more difficult calculations related to trigonometric and algebraic functional sets appearing in the computing process, including exp, ln, sin, and atan functions. In contrast, GMDH used polynomial quadratic equations. From this perspective as well, it seems the computational effort was less than GEP. Finally, the comparison results of ELM with MARS indicated that both algorithms had almost identical results.


## 4. Summary and Conclusions

Stilling basins are used widely as hydraulic dissipation structures in large dams. An accurate estimation ofCP′ within stilling basins is a critical issue for hydraulic engineering for the design of stilling basins. This study employed three soft computing methods, including ELM, GMDH, and M5MT, to estimate the CP′ that occurred during the hydraulic jump in the sloping channels.Different soft computing models were developed to estimateCP′ according to the dimensionless parameters and experimental data. ELM showed the lowest error values of RMSE, MAE, and SI in the training and testing stages. In addition, ELM has the highest correlation coefficient and Ia values than those obtained from GMDH and ELM. Therefore, the proposed soft computing models provided sufficiently accurate results due to this problem’s complexity.ELM formulated the matrices of coefficients for the prediction ofCP′. In addition, GMDH provided mathematical quadratic polynomials and combined input variables for the prediction of CP′. M5MT generated four simple rules for the estimation of CP′. It seems the application of M5MT was the most straightforward method for the prediction of CP′. It is noteworthy that complexity degree equations were derived from ELM and GMDH for the prediction CP′ have been more than M5MT rules.In summary, the ELM method provided a weight matrix for predictingCP′. GMDH generated second-order polynomial equations for the prediction of CP′. M5MT developed piecewise multiple linear regression equations for the calculation of CP′. GMDH and M5MT methods generated explicit and clear mathematical expressions to estimate CP′. In addition, M5MT divided the problem domain into subdomains and fitted local linear models to compute CP′. It seems that in terms of the degree of complexity of the developed models in estimating CP′, the ELM model has the highest complexity, followed by GMDH and M5MT models. However, in contrast to the complexity degree of the proposed models, their computational accuracy has increased for the estimation of CP′ so that the ELM model has the highest accuracy, followed by GMDH and M5MT.Finally, the comparisons of the results with previous research revealed that the proposed applications of soft computing methods have good performance for predictionCP′.Further works can be considered to use pressure field data for modeling with soft computing methods. It is also recommended that hybrid data-driven models with evolutionary algorithms be used instead of stand-alone data-driven models to figure out the coefficient of pressure fluctuations.



---

*Source: 2495631-2022-09-23.xml*