Hierarchical Fuzzy Systems Integrated with Particle Swarm Optimization for Daily Reference Evapotranspiration Prediction: a Novel Approach

Reference evapotranspiration (ET0) is a crucial element for deriving irrigation scheduling of major crops. Thus, precise projection of ET0 is essential for better management of scarce water resources in many parts of the globe. This study evaluates the potential of a Hierarchical Fuzzy System (HFS) optimized by Particle Swarm Optimization (PSO) algorithm (PSO-HFS) to predict daily ET0. The meteorological variables and estimated ET0 (using FAO-56 Penman–Monteith equation) were employed as inputs and outputs, respectively, for the PSO-HFS model. The prediction accuracy of PSO-HFS was compared with that of a Fuzzy Inference System (FIS), M5 Model Tree (M5Tree), and a Regression Tree (RT) model. Ranking of the models was performed using the concept of Shannon’s Entropy that accounts for a set of performance evaluation indices. Results revealed that the PSO-HFS model performed better (with Entropy weight = 0.93) than the benchmark models (Entropy weights of 0.77, 0.74, and 0.90 for the FIS, RT, and M5Tree, respectively). Furthermore, the generalization capabilities of the proposed models were evaluated using the dataset from a test station. Generalization performances revealed that the models performed equally well with the unseen test dataset and that the PSO-HFS model provided superior performance (with R = 0.93, RMSE = 0.59 mm d−1 and IOA = 0.94) while the RT model (with R = 0.82, RMSE = 0.90 mm d−1, and IOA = 0.83) exhibited the worst performance for the test dataset. The overall results imply that the PSO-HFS model could effectively be utilized to model ET0 quite efficiently and accurately.


Introduction
Irrigating crops to enhance agricultural productivity essentially require sufficiently large volumes of freshwater. Water-saving through carefully managed irrigation practices can be achieved by precise quantification of evapotranspiration (ET), which is used to develop correct irrigation scheduling, determine hydrologic water balances, simulate crop yields, and allocate water resources (Kisi 2016). Being an essential component of water balance, ET plays a vital function in controlling interactions among the atmosphere, soil, and vegetation . The measurement of ET may be performed through experimental methods including the Bowen ratio energy balance approach, eddy-covariance techniques, and lysimeter methods (direct methods) (Martí et al. 2015). Alternatively, ET can be obtained indirectly by calculating reference or potential evapotranspiration (ET 0 ) using climatological data. This indirect method has become popular in many parts of the world where direct measurements are not available or affordable due to complexity or costliness (Allen et al. 1998). As a universal approach of ET 0 estimation, the FAO-56 Penman-Monteith (FAO56-PM) (Allen et al. 1998) method has been recognized as the widespread reference method that can be employed in areas with varying ecological and climatic circumstances with no requirement of regional adjustments (Allen et al. 1998). FAO56-PM equation can be utilized to estimate ET 0 , which together with crop coefficient value provides an estimate of ET for a particular crop. Fewer recent studies also attempted to provide reliable estimates of ET 0 using a simplified version of the FAO56-PM equation (e.g., in Alavi and Rahimikhoob 2016;Zarei et al. 2021). The surface energy balance algorithm for land (SEBAL) approach using remote sensing data for ET 0 estimation has also been proposed in recent literature (Elkatoury et al. 2020).
In recent years, artificial intelligence or machine learning-based models have effectively been employed to model ET 0 in different hydrogeologic conditions. These models can map the complex and nonlinear relations between the input and output data quite effectively and accurately. Various models have been used in ET 0 modelling; among them, Artificial Neural Network (ANN) (Kumar et al. 2002; Gocić and Amiri 2021) models were the first implementation of machine learning tools to estimate ET 0 . Other recent implementation of machine learning tools in ET 0 modelling includes the usage of Adaptive Neuro Fuzzy Inference System (ANFIS) (Kisi and Zounemat-Kermani 2014;Petković et al. 2020;, Gaussian Process Regression (GPR) (Karbasi 2018), Wavelet-GPR (Karbasi 2018), Gene-Expression Programming (GEP) (Wang et al. 2019), M5 Model Tree (M5Tree) (Kisi 2016), Multivariate Adaptive Regression Splines (MARS) (Kisi 2016), Random Forest (RF) (Wang et al. 2019;Ferreira and da Cunha 2020;Salam and Islam 2020), Long Short-Term Memory Networks (Roy 2021), and Support Vector Machine (SVM) (Ferreira and da Cunha 2019; Chia et al. 2020;Salam and Islam 2020). Generally, artificial intelligence-based models provide superior performances over the typical equations in estimating ET 0 with similar datasets (Reis et al. 2019).
Although several machine learning-based ET 0 models were proposed recently, most of the modelling tasks (i.e., training, validation, and testing) were performed on the same station's meteorological data. There are few studies evaluating the model's performance with data gathered from other stations (i.e., from outside the training station). As an example, Wang et al. (2019) evaluated the potential generalization ability of RF and GEP models to estimate ET 0 using the data from several weather stations in China. The present research proposes an approach in which model training, validation, and testing are conducted in a single station, and the developed models are tested for the meteorological data from a nearby station to demonstrate the generalization ability of the proposed models.
Among various machine learning algorithms, tree-based algorithms such as RFs, Regression Trees (RT), and M5Tree have recently been gained significant attention due to their simplicity, robustness, and capability to provide accurate predictions of ET 0 (Chen et al. 2020). On the other hand, machine learning models derived from the theory of fuzzy logic have recently been utilized as an effective prediction system in various water resources management issues (Kord and Moghaddam 2014). Although an ANFIS, a variant of Fuzzy Inference System (FIS) has been successfully applied in developing ET 0 prediction models, the usage of ANFIS is hindered by the computational burden arising from a large number of rule bases especially for problems with larger input variables. This limitation is attributable to the fact that the number of rules in a fuzzy system escalates exponentially with the number of variables inputted into the system. Larger rule bases make the learning and fine-tuning of the rules and membership function parameters extremely challenging. In addition, larger rule bases reduce the generalization capability of tuned fuzzy systems when the number of training data is a limiting factor. To overcome this issue, FISs may be represented as a tree of smaller interrelated and interconnected FIS objects known as Hierarchical Fuzzy Systems (HFS), where the predictions from the lower-level FISs are utilized as predictors to the higher-level FISs making the fuzzy tree-based HFS computationally more efficient than a single monolithic FIS object.
HFS, an improved version of decision trees, provides reliable modelling using the concept of fuzzy logic principle. Although applied quite successfully in various research domains (Zheng et al. 2019), HFS models received extremely little attention in the fields of hydrological and agricultural research. Few recent studies related to hydrology and water resources management focused on the use of fuzzy logic-based decision trees (FDT). For instance, Sikorska-Senoner and Seibert (2020) employed an FDT instead of traditional trend analysis for quantifying the magnitudes and frequencies in the time series of floods. Wei and Hsu (2008) performed a comparison between three types of decision trees: neural decision trees, conventional decision trees, and FDTs to derive operating rules for a reservoir operation system. Their comparison results demonstrated the superiority of FDTs over the other two types of decision trees. Han et al. (2002) addressed uncertainty in real-time flood forecasting using FDTs. They concluded that although FDTs did not perform as good as the ANN models for river flow modelling in the test case, the glass box nature of fuzzy tree modelling could allow several valuable insights into the hydrological processes.
Although used extensively in various research domains, the classic machine learning-based models often suffer from overfitting and computational efficiency issues (Kisi et al. 2019). These drawbacks can effectively be solved through parameter tuning of machine learning models using evolutionary algorithms. For instance, Azad et al. (2018a) reported the potential applicability of Genetic Algorithm (GA), Particle Swarm Optimization (PSO), and Differential Evolution (DE) for tuning ANFIS parameters to provide estimations of water quality parameters. Azad et al. (2018b) utilized PSO tuned ANFIS for predicting multistep ahead (1-day to 5-day) river water flow. Kisi et al. (2017) reported that GA and PSO-based hybridized ANN models could model water level fluctuations quite effectively. In a different study, Azad et al. (2019a) proposed a hybridized GA-ANFIS model to predict air quality. Azad et al. (2019b) used four optimization algorithms, e.g., DE, GA, PSO, and Ant Colony Optimization (ACO) for obtaining optimal ANFIS parameters to model monthly precipitation. Later, Azad et al. (2021) assessed the abilities of ACO and GA in finding optimal parameters of ANFIS and Least Square SVM for monthly rainfall predictions.
Apart from that, several bio-inspired optimization tuned machine learning algorithms were believed to have an improved prediction accuracy related to ET 0 prediction. For instance, Alizamir et al. (2020) modelled ET 0 using GA and PSO optimized ANFIS models and compared their accuracy with the traditional ANFIS and ANN models.  utilized Biogeography-based Optimization (BBO), Firefly Algorithm (FA), PSO, and Teaching Learning Based Optimization (TLBO) to tune the ANFIS parameters for predicting ET 0 in subtropical climatic zones. SVR models optimized with Whale Optimization Algorithm (WOA) were found superior to the traditional SVR models in predicting ET 0 in three different climatic regions of Iran (Mohammadi and Mehdizadeh 2020). In another study, Tikhamarine et al. (2020a) proposed hybridized SVR models coupled with WOA, Multi-verse Optimizer (MVO), and Ant Lion Optimizer (ALO) for estimating monthly ET 0 at two meteorological stations in the north of Algeria with three different input scenarios. Grey Wolf Optimizer (GWO) tuned SVR model was also proposed by Tikhamarine et al. (2020b) for estimating monthly ET 0 at three meteorological stations in Algeria. Wu et al. (2020) coupled ELM with WOA and Flower Pollination Algorithm (FPA) to develop hybrid ELM models to estimate monthly pan evaporation and concluded that FPA significantly enhanced the accuracy of the ELM model. In a study conducted in an arid area in China (Zhu et al. 2020), the PSO tuned ELM model was found superior to traditional ELM, ANN, and RF-based models for predicting daily ET 0 using a limited number of meteorological information as inputs.
As far as the recent literature is concerned, HFS models have not yet been used in hydrological and agricultural research, especially in ET 0 modelling. Besides, the generalization capability of the developed ET 0 model for the target stations has not been evaluated using the HFS-based ET 0 modelling approach. A generalized prediction model which can be spatially validated by data from different stations eliminates the need to construct and calibrate a model specific to each station. Therefore, a general model that reflects a hydrological setting can be developed. Moreover, several bio-inspired optimization algorithms tuned soft computing approaches were believed to have an improved prediction accuracy related to ET 0 prediction (Tao et al. 2018;Wu et al. 2020). Therefore, it is plausible that optimization algorithms can improve the model performance significantly; therefore, in the current study, an optimization algorithm (e.g., PSO) was used to tune the modifiable parameters of individual fuzzy systems within the proposed HFS model (PSO-HFS).
Considering the importance of reliable estimates of ET 0 , the purposes of this study are to (1) assess the potentiality of PSO tuned HFS model (PSO-HFS) to predict daily ET 0 ; (2) weigh against the prediction capability of the proposed PSO-HFS with that of two tree-based models (RT and M5Tree) and a fuzzy logic-based model, FIS; (3) rank the proposed models with respect to their prediction accuracies utilizing several performance evaluation indices; and (4) evaluate the generalization capability of the proposed models outside the training station using data from a test station. According to the authors' understanding, this study is the first effort an evolutionary algorithm-tuned fuzzy decision tree (PSO-HFS) is employed to provide daily ET 0 prediction.

Materials and Methods
The study proposed a PSO-HFS model for predicting daily values of ET 0 from the inputoutput relationships of meteorological variables and ET 0 . Prediction of the proposed PSO-HFS model was then compared with that of three machine learning algorithms: a fuzzy logic-based FIS model and two tree-based models (RT and M5Tree). Comparison of prediction performances was evaluated using several statistical indices within the framework of Shannon's entropy that incorporated three benefits (higher values indicate better model performance: Correlation Coefficient, Nash Sutcliffe Efficiency Coefficient, and Index of Agreement) and three cost indices (smaller values indicate better model performance: Root Mean Squared Error, Mean Absolute Error, and Median Absolute Deviation) in the decisionmaking process. The proposed methodology was evaluated using the daily climatological data acquired from a weather station located in the Gazipur Sadar in Bangladesh. The developed models were then validated using daily climatic data from Ishurdi meteorological station in Bangladesh. A brief explanation of methodology parts is presented in the subsequent subsections.

Study Area and the Dataset
Meteorological variables were acquired from two different weather stations situated in the Gazipur Sadar Upazila of the Gazipur district and Ishurdi Upazilla of the Pabna district in Bangladesh. The weather station in Gazipur Sadar is positioned between 24.00°N latitude and 90.43°S longitude with an altitude of 8.4 m above the mean sea level. Meteorological variables including solar radiation, relative humidity, minimum and maximum temperatures, and wind speed were obtained for 15.5 years (from 01 January 2004 to 30 June 2019). Descriptive statistics of the meteorological variables for the training station are given in Table 1. It is perceived from Table 1 that the climatological variables demonstrated left (negative) skewness, which indicates that the distribution of data for all variables had an extended left tail than the right tail. Kurtosis, on the other hand, had both positive and negative values indicating that the datasets had both "heavy-tailed" (positive values of kurtosis) and "light-tailed" (negative values of kurtosis) distributions.
The data for the test station was acquired from 01 June 2015 to 31 December 2020 (2021 daily records of meteorological data and calculated daily ET 0 ). Statistical performance indices were computed for the entire (2021 entries: from 01 June 2015 to 31 December 2020), first half (1021 entries: from 01 June 2015 to 17 March 2018), and the second half (1020 entries: from 18 March 2018 to 31 December 2020) of the dataset for the test station. The selection of three sets of data allows investigating a better generalization capability of the model. Descriptive statistics of the meteorological variables of the test station are presented in Table 2. The locations of the weather stations in the study areas are presented in Fig. 1.
Meteorological variables obtained from the study areas across the period of study were utilized to estimate daily ET 0 by employing the FAO56-PM equation. These computed daily values of ET 0 and the meteorological variables were employed as model outputs and inputs, respectively for the proposed PSO-HFS and other models. This indirect approach of ET 0 estimation from meteorological variables has been widely accepted in circumstances when ET 0 estimates are extremely hard to obtain directly (Allen et al. 1998). The FAO56-PM equation is represented by: where, ET 0 represents reference evapotranspiration, mm d −1 ; R n is the net radiation at the crop surface, MJ m −2 d −1 ; G is the heat flux density of soil, MJ m −2 d −1 ; Δ is the slope of the saturation vapor pressure curve, kP a • C −1 ; γ is the psychometric constant, kP a • C −1 ; e s is the saturation vapor pressure, kP a ; e a is the actual vapor pressure, kP a ; u 2 is the wind speed at a height of 2 m, m s −1 ; and T mean is the mean air temperature at 2.0 m height, °C.
For the training station (Gazipur Sadar), computed ET 0 values ranged between 0.92 and 8.02 mm d −1 with the mean and standard deviation values of 3.80 and 1.32 mm d −1 , respectively. The distribution of ET 0 time-series had an extended right tail compared to the left tail as indicated by a positive skewness estimate of 0.30. The negative kurtosis estimate of -0.67 indicates a "light-tailed" distribution for the computed ET 0 values of the train ststion. On the other hand, the mean, standard deviation, skewness, and kurtosis values of the computed ET 0 for the entire dataset of the test station were found to be 3.67 mm d −1 , 1.24 mm d −1 , 0.28, and -0.62, respectively. For the first half of the dataset, the values were 3.57 mm d −1 , 1.25 mm d −1 , 0.35, and -0.62, respectively. The second half of the dataset comprised the following values of ET 0 : mean = 3.76 mm d −1 , standard deviation = 1.23 mm d −1 , skewness = 0.22, kurtosis = -0.59.

Proposed ET 0 Prediction Model: Hierarchical Fuzzy Systems (HFS)
Fuzzy Inference Systems (FIS) are regarded as one of the most effective tools for modelling dynamic and nonlinear systems with single output and multiple inputs (Takagi and Sugeno 1985;Sugeno and Yasukawa 1993). In general, however, the computational efficiency during the FIS training largely relies on the quantity of inputs to the FIS system and the quantity of rule sets, which generally increase exponentially as the number of input variables increases. A significant amount of rule sets not only reduces the computational efficiency but also creates difficulty in the tuning process of the rule base and membership function parameters. Moreover, an increased number of rule bases reduces the generalization capability of tuned FISs, especially in situations when the amounts of training data are scarce as it may be common in many hydrological applications. As a solution to these problems, an HFS consisting of smaller interconnected FIS objects can be implemented instead of a single massive FIS object with many input variables. In an HFS, the FISs are organized in 'hierarchical tree structures' in which the predictions from the lowerlevel FISs are employed as predictors to the higher-level FISs. With a similar number of input variables, an HFS usually requires fewer computational efforts compared to a single FIS. Fuzzy tree structures, based on which an HFS is constructed, can be of three major types for many practical applications: (a) incremental, (b) aggregated, and (c) cascaded or combined that combines both incremental and aggregated structures (Siddique and Adeli 2013). As cascaded tree structures are better suited for applications with both uncorrelated and correlated input variables, a cascaded or combined fuzzy tree structure was utilized in this research for constructing an HFS.
The first step of developing an HFS involves the creation of several FIS objects using the available input variables ranked based on their correlations with the output variable (ET 0 ). Both positively and negatively correlated input attributes were used to incorporate both the positive and negative impacts of input attributes on the output (ET 0 ) for prediction. Next, in T mean +273 u 2 e s − e a Δ + 1 + 0.34u 2 the second step, the input attributes were paired according to their ranks to create individual FIS objects as follows: • fis1: Maximum Temperature and Relative Humidity • fis2: Wind Speed and Sunshine Duration • fis3: Minimum Temperature and Sunshine Duration Then, the HFS as shown in Fig. 2, was constructed using the principle of FIS tree structure (Mathworks 2021). The constructed HFS had five two-input and one-output FIS objects (fis1, fis2, fis3, fis4, and fis5 in Fig. 2) of which the first three FISs (fis1, fis2, and fis3) received the ranked input attributes directly and produced intermediate predictions of ET 0 . The intermediate ET 0 values were integrated utilizing the remaining two FISs (fis4 and fis5). The HFS models were developed in the MATLAB environment.

Fuzzy Inference System (FIS)
The underlying principles of FIS are derived from fuzzy set theory that has received substantial interest over the past few years. FISs have successfully been applied as a reliable computing framework in various research areas (Jang et al. 1997). FISs have been regarded as effective prediction tools for modelling nonlinear processes due to their capability of accurately mapping the relationships (usually nonlinear) between input and output variables (Takagi and Sugeno 1985;Sugeno and Yasukawa 1993). A Sugeno Sugeno-type FIS (Sugeno 1985), also known as a Takagi-Sugeno-Kang FIS, is particularly better fitted for nonlinear system modelling. A Sugeno FIS builds input-output relationships through interpolating the outputs from multiple linear models. The basic structure of a Sugeno FIS consists of a rule base, a database, and a reasoning mechanism. The working principle of a basic FIS with three inputs, one output, and four rules is illustrated in a block diagram as shown in Fig. 3.
Rule bases of a FIS consist of fuzzy if-then rules, the database determines the types and numbers of MFs utilized in fuzzy rules, and finally, the reasoning mechanism accomplishes the fuzzy inference process (Jang et al. 1997). Several fuzzy if-then rules are utilized in a fuzzy inference process for producing a nonlinear mapping of input and output variables. A fuzzy rule consists of two parts: (a) the antecedent part of any rule specifies a fuzzy region within the input variable space, and (b) the consequent part specifies a fuzzy region within the output variable space. A Sugeno-type FIS introduced in 1985 (Sugeno 1985) was developed and utilized in this effort. The input and output MFs of the utilized Sugeno FIS were Gaussian and linear, respectively. The FIS-based models were developed in the MATLAB environment.

M5 Model Tree (M5Tree)
The development of M5Tree is based on the principles of the M5 method (Quinlan 1992). This method builds single trees corresponding to a method called 'M5', which makes use of the 'Standard Deviation Reduction' criterion. Model trees (MT) are combinations of traditional regression trees, which possess linear regression functions at the leaf nodes. Pruning and smoothing operators determine the contents of a leaf node for an MT, which is nothing but an RT without 'pruning' and 'smoothing' operations. MTs (Quinlan 1992) are machine learning algorithms that have demonstrated their predictive capabilities in various research domains (Bhattacharya and Solomatine 2005). MTs are 'inverted trees', i.e., root nodes are situated at the upper portion of the tree whereas the bottom portion of the tree contains the leaves.
MTs and RTs are the variants of Decision Trees (DT) that have been established for solving regression tasks (Quinlan 1992). However, MTs differ from RTs: while MTs generate linear models at their leaves, the RTs yield a constant value at their leaves. The linear models developed at the leaves are used to contain input-output relationships, which are then utilized to predict outputs for a given set of data. MTs are more efficient than RTs in handling large datasets and producing more accurate predictions. M5Tree utilizes 'divide-and-conquer' method that allows dividing the entire data space into smaller data sub-spaces (Bhattacharya and Solomatine 2005). In this approach, the input parameter space is narrowed down to several subspaces each of which represents a linear regression model. This unique data splitting procedure enables M5Tree to produce a hierarchial model tree that contains 'splitting rules' in its 'non-terminal nodes' and has 'expert models' in its leaves.

Regression Tree (RT)
RTs are decision trees that build simple, flexible, and easily interpretable models developed using input-output training patterns. RTs are associated with the principle of the 'Classification and Regression Tree (CART)' algorithm (Breiman et al. 1984;Krzywinski and Altman 2017). The CART algorithm follows three major stepwise procedures in building models: (a) building a complex tree, (b) pruning, and (c) selecting an optimal subtree. In the first step, a complex full tree with several terminal nodes is built using a binary split procedure. The complex tree built in the first step is pruned in the second step to prevent or at least reduce model overfitting. In the third step, an optimal subtree is selected by the CART algorithm to ensure the quality of prediction for new samples. The developed RTs deliver a predicted response by following the decisions from the beginning node (root node) to the leaf node within the tree. The leaf node of an RT contains the responses or outputs. The RT-based ET 0 prediction models were developed in the MAT-LAB environment.

Input-Output Datasets for ET 0 Prediction Models
The meteorological variables and computed ET 0 formed input-output datasets for the proposed ET 0 prediction models. The dataset of the training station comprises 5660 daily records (from 01 January 2004 to 30 June 2019) of meteorological variables and estimated ET 0 values. The entire dataset of 5660 entries was divided into train, validation, and test datasets: first 80% of the total data (4528 records: from 01 January 2004 to 24 May 2016) was employed to train and validate the proposed models while the last 20% (1132 records: from 25 May 2016 to 30 June 2019) was employed to test the developed models. The first 80% of the sequential data were randomized to minimize the effect of trends during the model training and validation processes. The randomized data were then split into two equal-sized datasets for training (first 40%) and validation (remaining 40%) sets. It is noted that the test set (last 20% of the entire dataset) was kept in sequence as this dataset was used to test the developed models for the actual nature of data. This technique of data partitioning allows better performance assessment for the constructed models (Francone 2001). For performance assessment, numerous statistical indices were computed on the test dataset.

Parameter Tuning of HFS: Particle Swarm Optimization
A parameter tuning approach was adopted to achieve optimum performance of the constructed HFS. Both the rule bases and input-output Membership Functions (MF) were tuned in two subsequent phases. In the first phase, learning of the rule bases was accomplished while input and output MF parameters were kept constant. After learning the new rules through the first phase, the parameters of the input-output MF, as well as the learned rules, were tuned simultaneously in the second phase. To achieve computational efficiency in the parameter tuning process, the tuned rule base obtained from the first phase was used as the initial condition for the second phase. This allows a fast parameter tuning process and quick convergence to global optima. Particle Swarm Optimization (PSO) has recently acquired substantial attention for having several useful features including less complex structure, robustness, flexibility, and simple usage. These beneficial attributes of PSO have facilitated parameter tuning processes of various data-driven models. Therefore, this study utilized PSO in both phases of parameter tuning to obtain optimal parameter values of the constructed HFS.
The PSO (Kennedy and Eberhart 1995), a swarm-based stochastic search algorithm, is encouraged by communal and psychological tenets. PSO is one of the most reliable optimization algorithms for parameter tuning of intelligent models as it has high computational efficiency, rapid convergence capability, little probability of being trapped into local optima, and easy implementation. Despite several advantages, PSO may prove to be inefficient for problems with high complexity in which the probability of being trapped onto local optima may increase. The PSO is connected to swarm intelligence principles, which feigns the communal characteristics of predation commonly observed in fish schooling or bird flocking. PSO considers each particle as a possible solution inside the search domain of an optimization issue. In contrast, the flight attributes of particles are regarded as the exploration method for the entire individuals within a community. In PSO, the dynamic update of a particle's velocity is determined by the particle's previous optimal location and the swarm population.
PSO considers the values of the particle's goal function to be the corresponding fitness values. These fitness values are used to calculate the particles' optimal position. The fitness values are also utilized to update the particles' past most advantageous location and the swarm population's optimum location. Thus, the PSO algorithm's control parameters determine the convergence of particles trajectories (Kennedy and Eberhart 1995). The PSO method converges by keeping records of each particle's best fitness values, finding the global best particle, and updating the locations and velocities of each particle. In the event that the convergence is not achieved, the iterative process continues until the optimization problem converges to its optimal solution or until the user-defined maximum number of iterations is satisfied. The PSO tuned HFS, proposed in this study, is referred to as a PSO-HFS.

Implementation of the Proposed PSO-HFS to Model ET 0
In developing the proposed models, ET 0 was considered as the target variable whereas the meteorological variables were used as the inputs. The inputs to and outputs from the PSO-HFS can be represented in the generalized form as: Determining the optimal parameter values, which have a significant impact on the output variable, is the most important step in developing predictive models. Parameter tuning was performed using a swarm-based optimization algorithm, PSO described earlier in the previous section (Sect. 3: Parameter tuning of HFS: Particle swarm optimization). To further improve the HFS model's accuracy, all tunable parameters (rule bases and parameters of both input and output MFs) were optimized using the PSO algorithm. To obtain optimal parameter values for both phases of a simulation, PSO was used with a maximum number of 500 iterations. The parameter tuning process was accomplished before reaching this maximum number of iterations. The optimal option sets of the PSO algorithm used in this (2) Reference evapotranspiration, ET 0 = f (meteorological variables) study are as follows: function tolerance = 1 × 10 -6 , inertia range = [0.1, 1.1], initial swarm span = 2000, minimum adaptive neighborhood size = 0.25, self-adjustment weight = 1.49, social adjustment weight = 1.49, and swarm size = 100. These optimal values were obtained upon performing several trials concerning a tradeoff between prediction precision and computational efficiency.
Eventually, the input-output datasets used to construct the prediction models were split into train, validation, and test sets. Additionally, FIS, M5Tree, and RT-based ET 0 prediction models were also developed solely for comparison purposes.

Ranking of Models: Shannon's Entropy
Shannon's entropy (Shannon 1948) was applied to assign weights to individual ET 0 prediction models which were ranked according to the weights assigned to them. Shannon's entropy incorporates a set of benefit (the higher values indicate better model accuracy) and cost (the lower values indicate better model accuracy) performance evaluation indices. The detailed calculation steps of Shannon's entropy can be found in , and are not repeated here.

Performance Evaluation Criteria
The proposed modelling approach was evaluated with several performance indicators as follows: Root Mean Squared Error, RMSE: Normalized RMSE, NRMSE: Accuracy: Mean Bias Error (MBE): Nash-Sutcliffe Efficiency Coefficient (NS):

Results and Discussion
The values of FAO56-PM estimated ET 0 were considered as the benchmark for evaluating the performances of the proposed PSO-HFS, and other models (FIS, M5Tree, and RT) developed for comparison purposes. Ten statistical performance assessment indices were computed for both the calibration (Training and Validation datasets) and testing (applied dataset) phases of model building. Performance assessment indices were calculated on the FAO56-PM computed and model projected ET 0 values. Performances of the constructed models to the computed statistical indices for the calibration and testing phases are presented in the subsequent paragraphs.

Performance of the PSO-HFS During the Training Phase
The training phase of model building is regarded as the most important step in the development of any prediction model. To prevent model over-or underfitting, training performance is compared with the validation performance using the validation data. Training, as well as validation phases, were accomplished concurrently, and several evaluation indices were computed on the model projected and FAO56-PM estimated ET 0 values. Performances of the proposed PSO-HFS and other tree-based models during training and validation stages Table 3. It is evidenced from Table 3 that all performance indices showed a reasonably good performance of the proposed PSO-HFS model during training and validation stages as evidenced by the negligible difference in values of the performance evaluation indices between these two phases.
Although not as accurate as of the proposed PSO-HFS model, the FIS-based model performed equally well during the training and validation phases (Table 3). Training performances of RT and M5Tree were observed relatively better than their validation performances especially on the cost indices (MAD, MAE, MAPRE, RMSE, and NRMSE) as shown in Table 3. However, their performances on benefit indices (e.g., accuracy, IOA, NS, and R) were found to be almost similar which indicates reasonably decent performances on the benefit indices. Overall, the training and validation performances of the proposed PSO-HFS model were obtained superior to the other tree-based models. Nevertheless, although performed differently on different performance indices as well as on the training and validation phases, all models presented in Table 3 produced quite acceptable results in the context of prediction modelling.
Although RT and M5Tree suffered slightly from model overfitting as evidenced by the results obtained from the cost indices, the models produced satisfactory results concerning the benefit indices. This is also acceptable because prediction models often show contradictory performances, i.e., one model may be deemed suitable based on the RMSE criterion whereas the other model may be a good performer based on the R criterion (Müller and Piché 2011;Roy andDatta 2019, 2020). These conflicting characteristics of the constructed models necessitate the inclusion of several performance evaluation indices instead of a few indices utilizing a suitable decision theory for judging the performance of any prediction model. A decision theory incorporating Shannon's entropy-based weighting system was applied to judge the performance and to rank the developed models for the test performances of the trained and validated models (Subsection 7.3: Ranking of the developed ET 0 prediction models).

Performance of the PSO-HFS During the Testing Phase
The proposed PSO-HFS and other tree-based models were further tested using the test data which were employed neither for training nor for validation of the models, i.e., the models were tested using data from outside the training and validation datasets. The performances during the testing phase (with applied dataset) computed on the estimated FAO56-PM and model-predicted (calibrated and validated) ET 0 values using several performance evaluation indices are presented in Table 4. As the results indicate, the models performed equally well when compared to the training and validation phases. It is observed from Table 4 that although the performances were slightly poor during the testing phase when compared to the calibration and validation phases, the model performance during the testing phase was excellent in the context of prediction modelling. The testing results produced accuracy > 0.97, IOA > 0.99, NS > 0.99, R > 0.95, MAD < 0.2, MAE < 0.3, MAPRE < 8%, and NRMSE < 0.1 which indicate an excellent model performance. Models' performance is deemed excellent when the NS statistic value is greater than 0.8 (Gupta et al. 1999) suggesting an efficient performance of the developed models. Moreover, the NRMSE (or scatter index) values in the testing phase were 0.052 and 0.093, respectively for the best (PSO-HFS) and the worst (RT) model. These NRMSE values also illustrate the excellent performance of the developed models based on the criteria set in Heinemann et al. (2012) and in Li et al. (2013). According to them, model performance is said to be excellent when NRMSE value is lower than 0.1, good when NRMSE value is between 0.1 and 0.2, fair when NRMSE value is between 0.2 and 0.3, poor when NRMSE is greater than 0.3. In general, model performances for all the developed models were satisfactory as indicated by the lower values of MAD, MAE, MAPRE, and NRMSE together with higher values of accuracy, IOA, NS, and R. As reported in Table 4, the performance of the PSO-HFS was found superior for majority of the performance assessment criteria, except for IOA and MAD criteria (Table 4). Both PSO-FIS and M5Tree produced similar values of IOA (IOA = 0.999) whereas M5Tree performed slightly better than PSO-HFS when the MAD criterion was considered (MAD = 0.068 mm d −1 and 0.066 mm d −1 , respectively for PSO-HFS and M5Tree).
Based on the computed performance evaluation indices presented in Table 4, the proposed PSO-HFS model exhibited a high degree of precision. Although a direct comparison is not possible due to dissimilar study conditions (datasets and modelling approaches),  Alavi and Rahimikhoob (2016). The performance of the proposed PSO-HFS is also satisfactory based on the MAE criterion (MAE = 0.148 mm d −1 ), which is in reasonable accord with the recently published ET 0 modelling approaches (Walls et al. 2020;Sattari et al. 2020a, b). The IOA value (IOA = 0.999) presented in this effort is superior to the IOA value of 0.946 reported in Proias et al. (2020). The R, NS, IOA, MAD, RMSE, and NRMSE values produced by the PSO-HFS model in this study are well agreed with those presented in Roy et al. (2021) for ET 0 prediction using the FA-ANFIS model. The model performances were also evaluated based on scatter and regression plots presented in Figs. 4 and 5. ET 0 estimates of the four modelling approaches with the benchmark FAO56-PM method during the test period were illustrated in Fig. 4. As can be observed, the PSO-HFS predictions were nearer to the FAO56-PM estimated ET 0 values than the FIS and RT models. On the other hand, fluctuations of model predictions were remarkably close to each other in the case of the PSO-HFS and M5Tree. This confirms the calculated performance evaluation indices presented in Table 4.
Comparison of the regression plots obtained from the model predictions and FAO56-PM estimates for the test dataset is illustrated in Fig. 5. It is apparent from the regression plots that the PSO-HFS had fewer scattered predictions compared to the FIS and RT models, and was closely followed by the M5Tree. Regression plots confirmed the superior performance of the PSO-HFS model over the FIS, RT, and M5Tree.
Prediction accuracies of the proposed PSO-HFS and three other tree-based models presented as box plots of absolute errors (Fig. 6) reveal that the PSO-HFS model outperformed other prediction models and that the accuracy of prediction for the M5Tree was better than FIS and RT models. It is worth mentioning that the prediction accuracy of the M5Tree was equally well with the PSO-HFS model accuracy. However, the M5Tree had more high magnitude absolute error values when compared to the PSO-HFS model. Therefore, it is evidenced that the PSO-HFS model is considered as the superior performer, among others.
Prediction performances of the proposed PSO-HFS and other tree-based models were also evaluated using spider plots, also termed as 'radar' (Rankin et al. 2008), which can assess the performance of multifunctional systems. In a spider plot, relevant performance indices are chosen and assigned to an axis on a multidimensional plot. Displaying relevant data, a spider plot can be employed to evaluate the performance of any multifunctional entity, including performances of several prediction models on a particular performance evaluation index calculated using the actual and model-predicted data. In this effort, spider plots were drawn to illustrate performances of the developed prediction models on two benefits (Accuracy and R) and four cost indices (RMSE, MAE, MAD, MBE). The results are presented in Fig. 7 which manifest the superior performance of the PSO-HFS model over the M5Tree, FIS, and RT models.

Ranking of the Developed ET 0 Models
It is an obvious and well-established fact that prediction models behave differently in terms of prediction accuracies when different performance evaluation criteria are used to compute the prediction performances ). This contrasting behavior of prediction models needs to be resolved to provide unbiased suitability of an individual model. To resolve this contradictory behavior of models, Shannon's entropy-based decision theory was utilized to offer a performance ranking of the considered prediction models. This ranking approach made use of six performance indices, three of them were benefit indices (greater values designate . These benefit and cost indices were incorporated in calculating the weights for standalone prediction models. The entropy weights calculated using Shannon's entropy are presented in Fig. 8. It is apparent from Fig. 8 that Shannon's entropy-based decision theory determined that the PSO-HFS (entropy weight = 0.93) model had superior performance, followed by the M5Tree (entropy weight = 0.90), FIS (entropy weight = 0.77), and RT (entropy weight = 0.74) models. It can, therefore, be concluded that the PSO-HFS model achieved higher accuracy than the other benchmark models considered in this effort.

Generalization of Developed Models for a New Unseen Test Dataset
The PSO-HFS and other tree-based models developed at the training station (Gazipur Sadar) were validated using meteorological data obtained from a test station (Ishurdi station). Three distinct sets of data of the test station were inputted to the developed models for predicting daily ET 0 , which were then compared with the estimated ET 0, and different performance assessment indices were calculated using the model predicted and FAO56-PM estimated ET 0 values. The performance evaluation results in terms of various statistical indices are shown in Table 5. As the results indicate, the models performed equally well when compared to the results of the training station. The model performances were satisfactory based on the computed statistical indices: the model produced higher values of accuracy, NS, IOA, and R as well as lower values of RMSE, NRMSE, and MBE for all three datasets.

Conclusion
The potential of the PSO tuned HFS modelling approach (PSO-HFS) for the prediction of ET 0 using climatic variables was explored in this research. The study revealed that modelling of daily ET 0 can efficiently be performed using the Fuzzy logic-based PSO-HFS model, specifically when a Sugeno type FIS is employed for the construction of PSO-HFS. Five input attributes (climatic variables) such as solar radiation, relative humidity, minimum and maximum temperatures, and wind speed were utilized to predict the daily ET 0 . The PSO-HFS was constructed from five FIS objects built using the ranked input attributes (correlations of the input attributes with the output attribute, ET 0 ). The training and validation with the dataset of the train station revealed that the developed PSO-HFS adequately mapped the input-output patterns of the train station dataset. Therefore, a PSO-HFS can effectively be applied in predicting ET 0 with  for the test dataset. The proposed modelling tool provides a promising approach for ET 0 estimation in sub-tropical climates.
The study applied data from one weather station and the developed models were tested for the unseen test dataset obtained from a nearby station. It is worthwhile to assess the usability of the proposed PSO-HFS modelling approach by including weather stations with varying climatic zones. Future research may be directed towards exploring and comparing other bio-inspired optimization algorithms for the parameter tuning process of the HFS models.
Funding This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Availability of Data and Material
Datasets and other materials are available with the authors, and may be accessible at any time upon request.
Code Availability MATLAB codes are available with the first author.

Conflict of Interest
The authors declare that there is no conflict of interest.