Evaluation Method of Multiobjective Functions' Combination and Its Application in Hydrological Model Evaluation

Parameter optimization of a hydrological model is intrinsically a high dimensional, nonlinear, multivariable, combinatorial optimization problem which involves a set of different objectives. Currently, the assessment of optimization results for the hydrological model is usually made through calculations and comparisons of objective function values of simulated and observed variables. Thus, the proper selection of objective functions' combination for model parameter optimization has an important impact on the hydrological forecasting. There exist various objective functions, and how to analyze and evaluate the objective function combinations for selecting the optimal parameters has not been studied in depth. Therefore, to select the proper objective function combination which can balance the trade-off among various design objectives and achieve the overall best benefit, a simple and convenient framework for the comparison of the influence of different objective function combinations on the optimization results is urgently needed. In this paper, various objective functions related to parameters optimization of hydrological models were collected from the literature and constructed to nine combinations. Then, a selection and evaluation framework of objective functions is proposed for hydrological model parameter optimization, in which a multiobjective artificial bee colony algorithm named RMOABC is employed to optimize the hydrological model and obtain the Pareto optimal solutions. The parameter optimization problem of the Xinanjiang hydrological model was taken as the application case for long-term runoff prediction in the Heihe River basin. Finally, the technique for order preference by similarity to ideal solution (TOPSIS) based on the entropy theory is adapted to sort the Pareto optimal solutions to compare these combinations of objective functions and obtain the comprehensive optimal objective functions' combination. The experiments results demonstrate that the combination 2 of objective functions can provide more comprehensive and reliable dominant options (i.e., parameter sets) for practical hydrological forecasting in the study area. The entropy-based method has been proved that it is effective to analyze and evaluate the performance of different combinations of objective functions and can provide more comprehensive and impersonal decision support for hydrological forecasting.


Introduction
The parameters optimization of hydrological models has always been an important research content in the hydrological field, and it has a crucial impact on the overall performance of hydrological models and the reliable results of hydrological forecasting [1]. The essence of hydrological model parameters optimization is that adjusting the parameters of the hydrological model to address the errors and uncertainties in the watershed model simulation so that the model output values are as close as possible to the actual values [2,3].
The behavior and performance evaluation of the hydrological model is commonly made and reported through comparisons of objective function values of simulated and observed variables. Based on the comparison results, the hydrologist can make mathematical subjective and/or objective estimates of the "closeness" of simulated behavior of the model to observations made within the watershed [4].
The single objective function which can only consider one aspect of hydrological features is adopted in the traditional optimization methods for hydrological model parameter [5,6]. The engineering practices show that parameter optimization of the hydrological model is intrinsically a high dimensional, nonlinear, multivariable, combinatorial optimization problem which involves a set of different objectives (multiobjective) [7]. To fully excavate the various hydrological feature information contained in hydrological data, multiobjective research has become one of the important directions in hydrological model calibration [8].
Two or more conflicting objectives often have to be simultaneously optimized and the trade-off decisions should be taken between these objectives in the multiobjective optimization [9,10]. The multiobjective optimization model such as the Pareto theory [11] can reduce the degree of simplification of the model and can more fully reflect the requirements of engineering design, thus providing a decision basis for obtaining the optimal solution set.
Parameters optimization of hydrological model is a typical unconstrained multiobjective optimization problem which determines the model parameters' values by the optimization method to make the selected multiple objective functions achieve maximum or minimum simultaneously. Assuming that all objective functions need to be minimized, the multiobjective optimization problem of hydrological model parameters is expressed in the following equation [12]: f(X) � min f 1 (X), f 2 (X), . . . , f N (X) , where f i (X), i � (1, 2, . . . , N) is the i-th objective function; N is the number of the objective functions; X is a vector composed of decision variables, namely, the parameters combination of the hydrological model; and D is the number of hydrological model parameters to be optimized. The optimization result is the noninferior parameter combination sets in which feasible regions are limited by the values' ranges of the model parameters. Some researchers have applied the multiobjective method for the parameters optimization problem of hydrological models and achieved good results [13].
In studies and literature of hydrologic modeling, there are a large number of objective functions such as the Nash-Sutcliffe efficiency coefficient (NSE), root mean square error (RMSE), and coefficient of determination (R 2 ) that are frequently used. Because the objective functions are mutually constrained and each objective function places different emphasis on different systematic or dynamic behavioral errors, it is difficult for a hydrologist to clearly evaluate the advantages and disadvantages of the solutions of parameters optimization.
Most of the related research studies focus on the design of multiobjective algorithm and construction of objective functions, which provides important theoretical basis and practical experience for multiobjective optimization of hydrological model parameters. However, they do not provide much guidance for the selection of the actual objective function combination. Thus, how to study the influence of different objective function combinations on the optimization results of hydrological model parameters and make decisions on these Pareto solution sets has become one of the most important issues.
The technique for order preference by similarity to ideal solution (TOPSIS) which was first proposed by Hwang and Yoon is a well-known method for solving a multiple criteria decision-making problem. TOPSIS is designed based upon the concept that the chosen alternative should have the shortest distance from the positive ideal solution (PIS) and the farthest from the negative ideal solution (NIS) [14]. Thus, it can be adapted to evaluate and give the ranking order of all alternatives in solving a multiple attribute decision-making (MADM) problem. Because of its significant features of comprehensiveness, ease of programming, and reliability, the TOPSIS has been adopted in the decision-making of multiobjective optimization by computing the closeness of the different alternatives to an ideal solution [15].
In the state-of-the-art optimization algorithms, the artificial bee colony (ABC) algorithm [16] which was proposed by Karaboga in 2005 has the advantages of fast convergence speed, strong robustness, and so on. And it has been widely applied in different fields [17,18]. In this paper, various objective functions related to parameters optimization of hydrological models were collected from the literature and constructed to nine combinations. Then, a multiobjective artificial bee colony algorithm named RMOABC which adopts the mechanisms of the regulation operator and adaptive grid is employed to optimize the hydrological model and obtain the Pareto optimal solution [19]. The parameter optimization problem of the Xinanjiang hydrological model was taken as the application case for long-term runoff prediction in the Heihe River basin. Finally, the TOPSIS method based on the entropy theory is adapted to sort the Pareto optimal solutions to compare the combinations of objective functions and obtain the comprehensive optimal combination. The paper is organized as follows. The related literatures are reviewed in Section 2. Section 3 describes the hydrological Xinanjiang model, the study area, the RMOABC algorithm, and the TOPSIS method. Section 4 summarizes the objective functions and their combinations commonly used in the literature for hydrological model parameters optimization. The selection and evaluation framework of objective functions for hydrological model parameter optimization is proposed in Section 5. The experiment settings, results, and corresponding analyses are discussed in Section 6, and finally, the conclusions are drawn in Section 7.

Literature Review
In recent years, the ABC algorithm's advantages of great accuracy and satisfactory convergence speed make it suitable for solving the multiobjective optimization problems. Hedayatzadeh et al. designed a multiobjective artificial bee colony (MOABC) based on the Pareto theory and ε-domination notion [20]. The performance of a Pareto-based MOABC algorithm has been investigated by Akbari et al. on CEC'09 datasets [21], and the experimental results show that compared with the other multiobjective algorithms, the variants of multiobjective ABC can find solutions with competitive convergence and diversity within a shorter period of time.
The parameter calibration or optimization for the hydrological models has entered the era of multiobjective optimization, and a lot of literatures focus on the multiobjective research studies [22][23][24]. For the ABC algorithm, a novel multiobjective evolutionary algorithm named multiobjective artificial bee colony (MOABC) algorithm is presented and applied in long-term cascaded hydropower system dispatch in [25]. Pérez et al. presented a multiobjective artificial bee colony-based optimization approach to optimize the allocation of monitoring stations in the design of water quality monitoring networks in river basins [26].
In these studies, a number of objective functions were commonly used for a variety of model calibrations in hydrology such as mean squared error, absolute mean/maximum error, residual bias, Nash objective function [22], the sum of square errors [27], and so on. To evaluate and compare the objective functions, some scholars have also made some efforts. For instance, nine different efficiency criteria for the evaluation of model performance are described and compared through three examples involving an observed streamflow hydrograph from the Wilde Gera catchment in Germany [4]. And Guo et al. presented a method for the comparison of the influence of different objective function combinations on the optimization results of daily streamflow forecasting by a hydrological model [28].
In the evaluation and ranking methods, the TOPSIS is a well-known and widely used method for solving the multiple attribute decision-making problems. For example, Lia et al. extended the concept of TOPSIS to develop a methodology for solving a multiple objective decision-making problem [29]. Abo-Sinna et al. presented a method based on TOPSIS for selection of the large-scale multiobjective nonlinear programming problems with block angular structure [30]. Baky and Abo-Sinna extended the TOPSIS approach to solve bilevel MODM problems [31]. An iterative clustering around latent variable (CLV) based objective entropy weighted TOPSIS approach is proposed for benchmarking building energy performance in a multifactor manner [32]. And a multiobjective mathematical model under several constraints is developed to attain multiple goals simultaneously for sustainable energy investment planning [33].
The TOPSIS method is also adapted and applied for the multiobjective optimization problems (MOPs). For instance, Liu et al. proposed a new method of water distribution system design based on the multiobjective optimization algorithm, and the TOPSIS was used to rank the Pareto optimal solutions [34]. An intuitionistic fuzzy hybrid TOPSIS approach is proposed to handle risk assessment and is applied as a case study on a gas refinery for the welding and lamination task [35]. The concept of TOPSIS was extended to develop a methodology for solving multilevel nonlinear multiobjective decision-making (MLN-MODM) problems of maximization type [36]. TOPSIS was also used to select the optimal cutting conditions for trimming carbon fiber composite that could generate minimum tool wear, surface roughness, and tool temperature simultaneously while maintaining high production rates [15].
For the multiobjective parameter calibration of hydrological models, the TOPSIS has also been adopted by some scholars. A method which combines a genetic algorithm with TOPSIS for the Xinanjiang model is presented to handle the multicriteria parameter calibration problem with a single procedure [37]. A MOPSO algorithm and an entropy-based TOPSIS ranking method were employed to calibrate and verify the Xinanjiang model in Misai catchment [38].
In summary, as the hydrological model and the multiobjective optimization method are determined, the selection and combination of the objective functions become one of the most important factors affecting the prediction quality of the hydrological model. And the research studies on this field are still in the initial stage. Therefore, this paper mainly focuses on collecting the objective functions of hydrological models' optimization from the related literature, optimizing the parameters of the Xinanjiang hydrological model for river runoff forecasting in the Heihe River basin by the RMOABC algorithm and then evaluating and ranking the obtained nondominated solutions based on the TOPSIS method. It would provide a reference method and help for the selection and the combination of objective functions for parameters calibration of hydrological models.

Hydrological Model and Study Area.
In this paper, a modified Xinanjiang hydrological model [39] with two sources, surface runoff and underground runoff, was selected as the hydrological model. The Xinanjiang model is a conceptual rainfall-runoff (RR) model proposed by Renjun Zhao in 1973 and has been widely used in the humid and semihumid areas in China [40]. The model has ten parameters which have been proved to simulate the runoff well in the working practices [41], including the potential evapotranspiration ratio to pan evaporation K; the percentage of impervious in the catchment IMP; the exponential parameter with a single parabolic curve B, which represents the spatial distribution nonuniformity of the soil moisture storage capacity over the catchment; the average storage capacity of soil moisture of the upper layer WUM (unit: mm); the average storage capacity of soil moisture of the lower layer WLM (unit: mm); the average storage capacity of soil moisture of the deep layer (unit: mm); the coefficient of soil evaporation in the lower and deep layer of the river basin C; the steady infiltration rate FC (unit: mm/ hour); the outflow coefficients of the free water storage to groundwater relationships KKG; and the confluence parameters of surface runoff Kr. The value ranges of the ten parameters are shown in Table 1. The detailed description of the modified Xinanjiang model can be found in bib41 [7,12,41].
Computational Intelligence and Neuroscience 3 The Heihe River watershed in Northwest China was taken as the study area, and parameter optimization for runoff prediction which is a common concern in hydrology researches was taken as the case study. The observation data of daily rainfall, surface evaporation, and runoff from January 1,1990, to December 31, 1994, were selected for the experiments [7,12]. The period from January 1, 1990, to December 31, 1990, was chosen as the calibration period. The details of the data can be found in our previous papers [7,12,42].

The RMOABC Algorithm.
The RMOABC algorithm is proposed to extend the original ABC algorithm based on Pareto theory [43] for handling the multiobjective optimization problems [19]. The regulation operator mechanism is integrated to balance the local search and global search in the evolution process, and the adaptive grid [44] mechanism proposed in the PAES (Pareto archive evolutionary strategy) algorithm was utilized in the RMOABC to produce welldistributed nondominated Pareto solution set in the external archive. Each nondominated solution can be mapped in a certain location in the grid according to the values of multiobjective functions. The grid can adaptively maintain the distribution of candidate solutions stored in the external archive in a uniform way in the evolution process. The details of the RMOABC algorithm and how to integrate the algorithm with the Xinanjiang model can be found in the literature [12,19].

TOPSIS Method.
The technique for order preference by similarity to ideal solution (TOPSIS) was proposed by Hwang et al. in 1981, and it is a useful technique in dealing with multiattribute or multicriteria decision-making problems in the real world [22]. It is based on the concept that the chosen alternative should have the shortest distance to the positive ideal solution (PIS) and the furthest distance from the negative ideal solution (NIS). The positive ideal solution consists of the best performance values among any other alternatives for each attribute, whereas the negative ideal solution is the combination of the worst performance measures [45,46]. The TOPSIS method can produce a clear preference order of a set of competing designs, which can help decision makers carry out analysis, comparison, and ranking of the alternatives for the problems needed to be solved [22].
For the multiobjective optimization problems, TOPSIS can transfer the m-objective (criteria) problem, which are conflicting and difficult in comparison, into a two-objective (the shortest distance from the PIS and the longest distance from the NIS) problem. Then, the biobjective problem can be solved by calculating the relative closeness (overall performance coefficient) to the ideal solution [47]. Thus, TOPSIS can prioritize the nondominated solutions, rank the compromise solutions, and provide decision support for the researchers.

Multiobjective Functions for Calibration of Hydrological Model
The objective of model calibration is to select the model parameters so that the model simulates the hydrological behavior of the catchment as closely as possible. Objective functions which are defined as functions that must be minimized or maximized in the process of model parameter calibration are commonly used by hydrologists to provide a mathematical estimate of the "closeness" of the simulated behavior to the observed measurements [48]. Each objective function may place different emphasis on different types of simulated and observed behaviors; thus, the selection of the combination of the objective functions is crucial to the optimal parameter result and has been a challenge for the hydrologists.
In this section, we summarized the objective functions and their combinations commonly used in the literature for hydrological model parameters optimization. [49]. The Nash-Sutcliffe efficiency coefficient (NSE) [50] and its logarithmic form (LNNSE) [49] are the two common objective functions originating from statistical theory as the primary objectives for model calibration. The two objective functions are defined in equations (2) and (3), respectively. The NSE and the LNNSE functions are a pair of conflicting objective functions [7,12], where the NSE function evaluates the ability to reproduce all stream flows, but it is known as biased to predict peak flows. On the contrary, the LNNSE function emphasizes low flows.

Objective Functions in Literature
where N is the total number of time steps in the calibration period, Q sim,t is the simulated runoff data value at time t, Q obs,t is the observed data value at time t, and Q obs is the average of the observation data. The two objective functions are optimized to the maximum at the same time in the parameters optimization process. The larger value means the better performance of the model forecasting.  [51]. Madsen applied the following numerical performance statistics to measure the different calibration objectives: (1) Overall volume error: a good agreement between the average simulated and observed catchment runoff volume (i.e. a good water balance).
(2) Overall root mean square error (RMSE): a good overall agreement of the shape of hydrograph.
(3) Average RMSE of peak flow events: a good agreement of the peak flows with respect to timing, rate, and volume.
(4) Average RMSE of low flow events: a good agreement for low flows.
In Equations (4)- (7), N is the total number of time steps in the calibration period; Q obs,t is the observed discharge data value at time t; Q sim,t is the simulated runoff data value at time t; M p is the number of peak flow events; M l is the number of low flow events; n j is the number of time steps in peak flow event j; and n k is the number of time steps in low flow event k. Peak flow events are defined as periods where the observed discharge is above a given threshold level. Similarly, low flow events are defined as periods where the observed discharge is below a given threshold. The peak flow events were defined as periods with flow above a threshold value of 0.75 m 3 /s, and low flow events were defined as periods with flow below a threshold value of 0.2 m 3 /s. [52]. In literature [52], the mean squared logarithmic error E MSL and the mean fourth-power error E M4 of measured flow and forecasted flow were selected as the objective functions. The two objective functions are defined in equations (8) and (9):

Objective Functions in Literature
where N is the total number of time steps in the calibration period, Q sim,t is the simulated runoff data value at time t, and Q obs,t is the observed data value at time t. Due to the logarithmic operation, E MSL makes the fitting error of low flow more significant, so it is more focused on the fitting of low flow value. For the impact of the fourth power, the fitting error of high flow value has a greater contribution; thus, E M4 is more focused on the fit of high flow values. The experimental results showed that there is a typical noninferior relationship between the objective functions E MSL and E M4 . [53]. The parameters of the ANN-based monthly streamflow forecasting model are estimated by several multiobjective optimization algorithms with three pairs of objective functions. (1) mean squared logarithmic error (E MSL ) vs. mean squared derivative error (E MSD ), (2) mean fourth-power error (E M4 ) vs. mean squared derivative error (E MSD ), and (3) mean squared logarithmic error (E MSL ) vs. mean fourth-power error (E M4 ). The definition of E MSD objective function is shown in the following equation:

Objective Functions in Literature
where N is the total number of time steps in the calibration period, Q sim,t is the simulated runoff data value at time t, and Q obs,t is the observed data value at time t.
The E MSL function is more suitable for low flows due to the logarithmic transformation. E M4 is considered to be an indicator of the goodness of fit for high flows because it amplifies the effect of the deviation. The experiments demonstrated that E MSD can be taken as an indicator of the fit of the shape of the hydrograph, and it should be used in combination with E MSL or E M4 . [3]. The Nash-Sutcliffe efficiency coefficient (NSE) and the overall runoff relative error R E were selected as the objective functions in [3]. R E is defined in the following equation:

Objective Functions in Literature
where N is the total number of time steps in the calibration period, Q sim,t is the simulated runoff data value at time t, and Q obs,t is the observed data value at time t. [54]. There are four metrics used in [54]-two common statistical metrics that Computational Intelligence and Neuroscience emphasize high and low flows and two hydrological metrics that emphasize the water balance and the midrange flow regime. These metrics capture four different but important components of the hydrograph while allowing for a comparison between statistical and hydrological metrics. The first statistical metric is the commonly used root mean squared error (E RMS ) which has been defined in equation (5), and it emphasizes fitting the high flow portions of the hydrograph.

Objective Functions in Literature
The second statistical metric is transformed root mean squared error E TRMS . The simulated and observed flow time series are first transformed by a Box-Cox transformation (equation (12) with a λ value of 0.3, which has a similar effect as a log transformation. E TRMS emphasizes low flow portions of the hydrograph and is defined in equation (13).
where N is the total number of time steps in the calibration period, Z sim,t is the transformed simulated runoff data value at time t, and Z obs,t is the observed data value at the time. After the Box-Cox transformation, the observed and simulated flow data can help the parameter estimation method to better adapt to the low flows. The third metric is the runoff coefficient error (E ROC ) which captures the overall accuracy of the water balance by first combining the flows into one characteristic hydrological descriptor. It is defined in the following equation: where Q sim,annual and Q obs,annual are the simulated and observed mean annual runoff volume and P obs,annual is the mean annual precipitation which is 388.95 mm/year in the research area. The final metric is the slope of the flow duration curve error (E SFDC ), which measures how well the model captures the distribution of midlevel flows. The slope of a watershed's flow duration curve indicates the variability, or flashiness, of its flow magnitudes. E SFDC metric is calculated by equation (15), and it is simply the absolute error in the slope of the flow duration curve between the 30 and 70 percentile flows.
where Q sim, 30 and Q sim,70 are the 30 and 70 percentile flows of the simulated flow duration curve and Q obs, 30 and Q obs70 are the 30 and 70 percentile flows of the observed flow duration curve. [55]. The performance measures used in [55] are the Nash-Sutcliffe efficiency coefficient (NSE), the coefficient of determination for linear regression (R 2 ), and percent bias (P BIAS ). The coefficient of determination for linear regression is a measure of how well the regression line represents the data. The coefficient of determination R 2 is defined as the squared value of the coefficient of correlation which is defined in the following equation:

Objective Functions in Literature
where N is the total number of time steps in the calibration period, Q sim,t is the simulated runoff data value at time t, Q obs,t is the observed data value at time t, Q sim is the average of the simulated data, and Q obs is the average of the observation data. The range of R 2 lies between 0 and 1 which describes how much of the observed values can be reproduced by the prediction. A bigger value means better prediction. P BIAS which is defined in equation (17) measures the average tendency (deviation) of the simulated data to their observed counterparts and is expressed as a percentage. The optimal value of P BIAS is 0.0, with low-magnitude values indicating accurate model simulation [56]. [8]. Three evaluation functions are selected as the objective functions for the runoff prediction of the hydrological model: the root mean square error E MSR which represents the overall simulation result, the water volume balance function E RMOV , and the mean squared logarithmic error E MSL (overall volume error). E MSR is defined in the following equation:

Objective Functions in Literature
4.9. Objective Functions in Literature [28]. As shown in equations (19)-(24), six objective functions were constructed based on the absolute error of the model forecasted flow and the measured flow to establish three objective function combinations, and a simple method for the comparison of the influence of different objective function combinations was presented in [28]. The objective function combinations include E RMS L vs. E RMS H , R 2 L vs. R 2 H , and E MA L vs. E MA H . Among them, the subscripts H and L denote the flood peak flow objective and the low flow objective, respectively. 6 Computational Intelligence and Neuroscience where N is the total number of time steps in the calibration period; Q obs,t is the observed discharge data value at time t; Q sim,t is the simulated runoff data value at time t; Q obs is the average of the observation data; M p is the number of peak flow events; M l is the number of low flow events; n j is the number of time steps in peak flow event j; and n k is the number of time steps in low flow event k. Peak flow events or low flow events are defined as periods where the observed discharge is above or below a given threshold level, respectively. Similar to Section 4.2, the low flow threshold is 0.2 m 3 /s and the peak flow threshold is 0.75 m 3 /s. In the literature [28], the researchers found that the combination of E RMS L2 , E RMS H2 , and E MA L can effectively balance the multiple objectives and can provide more comprehensive decision-making support for the forecasters.

Selection and Evaluation Framework of Objective Function Combination
Solving a multiobjective optimization problem of parameter calibration for the hydrological model is to search for the nondominated solutions that can be expressed in terms of Pareto theory in the objective space of parameters. However, for the hydrologist, the nondominated solutions should be evaluated and ranked, and a set of parameters (i.e., one of the nondominated solutions) should be selected from the Pareto set for the application of the hydrological model. Thus, after the multiobjective optimization algorithm calibrating the parameters of a hydrological model and obtaining a set of nondominated solutions, the task of how to evaluate and rank these solutions corresponding to the objective function combination and select the optimal solution that meets the needs of decision makers has become more important. In this paper, a widely used ranking method called technique for order preference by similarity to ideal solution (TOPSIS) was used to prioritize the nondominated solutions and help the hydrologists find the appropriate combination of objective functions.

Entropy-Based TOPSIS.
TOPSIS approach was originally proposed for the comprehensive evaluation of multiple scenarios in the case of multiple attributes. It is one of the best-known ranking methods of alternatives for MADM problems and can be easily adapted to different fields. If the evaluation alternative is closest to the ideal solution and far away from the negative ideal solution, it is the best alternative; otherwise, it is the worst [46]. For the multiobjective parameters optimization problem of the hydrological model in this paper, each optimization objective function represents an attribute. Thus, the nondominated solutions generated from RMOABC algorithm were treated as alternatives and the objectives were taken as attributes.
In this paper, a TOPSIS method is employed to rank the nondominated solutions, and the entropy theory of information is also adapted to assign objective weights reasonably for different attributes. The entropy-based TOPSIS procedure consists of the following steps [32,38]: The decision matrix D is constructed according to the Pareto optimal solution set. Assume that there are n Pareto optimal solutions (alternatives) A i (i � 1, 2, . . ., n) to be evaluated against m criteria (objective functions) C j (j � 1, 2, . . ., m). Then, as shown in equation (25), the decision matrix D contains performance ratings for the i-th alternative (Pareto solution) A i with respect to the j-th criteria (objective function) C j , which are denoted as x ij .
In general, the evaluation criteria are classified into two types: benefit and cost. For criteria of benefit, larger value means more valuable, while for the criteria of cost criteria, smaller value means more valuable.
Then, the above decision matrix D is normalized and trended to obtain a normalized matrix R. According to the criteria type of benefit or cost, the normalized value r ij can be calculated as follows: Computational Intelligence and Neuroscience x ij are the maximum and the minimum values of the j-th criteria among n alternatives and r ij represents the normalized performance of A i with respect to attribute C j . After the normalized transformation, the normalized decisionmaking matrix is represented as follows: where 0 ≤ r ij ≤ 1 which means the optimal value for each attribute is 1 and the worst value is 0. This step can transform various attribute dimensions into nondimensional attributes to allow comparisons across criteria. (3) Calculate the weighted normalized decision matrix.
Considering the different importance of each criterion, the weighted normalized decision matrix Z can be constructed as follows: where w j is the weight of the j-th criterion, w j > 0, and m j�1 w j � 1. Generally, the weights of attributes are determined from the experience of decision maker depending on the specific problem. For more accuracy, a reasonable method based on the theory of information entropy for assigning objective weights of attributes is adapted to calculate objective weights [38]. According to the definition of entropy, the entropy value of each attribute (objective) j can be measured as follows: where H j is the entropy value of objective j, k is a constant (k � 1/ln(n)), n is the amount of nondominated solutions, and y ij is defined as where x ij is the i-th nondominated solution (alternative) under the j-th criteria (objective function). Thus, if the amount of objectives is m, the weight for each criterion j is given by the following equation: (4) Determine the positive ideal solution (best) and negative ideal solution (worst). The positive ideal solution (PIS), whose index consists of the best value for each criterion, is defined in equation (32). The negative ideal solution (NIS), whose index consists of the worst value for each criterion, is defined in equation (33).
(5) Calculate the separation measures by the m-dimensional Euclidean distance for each alternative. The separation from the positive ideal alternative can be expressed as equation (34) and similarly, the separation from the negative ideal alternative can be expressed as equation (35).
(6) Calculate the relative closeness (overall performance coefficient) to the ideal solution.
Finally, by sorting the closeness coefficients C * i calculated in Step 6, the ranking order of all alternatives can be obtained. According to the concept of the TOPSIS method, the higher value of C * i means that it has higher priority, and on the contrary, the smaller value of C * i means it has lower priority. Thus, the alternative with the highest C * i value will be the optimal alternative.

Combinations of Objective Functions.
In this section,we summarized and collated the objective functions and their combinations from the relevant literatures on hydrological model parameter optimization. Information on the combination of objective functions, related descriptions, objective function characteristics, and related literature is summarized in Table 2. The hydrological model parameter optimization experiments will be carried out according to these combinations of objective functions, and the runoff 8 Computational Intelligence and Neuroscience prediction of the hydrological model will be carried out. We have summarized a total of 19 objective functions and 9 objective function combinations whose number has been denoted from 1 to 9.
There are several reasons why we do not take all the object functions as the combination. (1) All objective functions are researched and summarized by many scientists in hydrological fields based on domain knowledge and practical work experience. It cannot be replaced and modified at will. (2) Because there are 19 kinds of objective functions, multiobjective functions cannot be optimized by algorithms at the same time, and an objective solution set that can balance each objective function cannot be obtained.
(3) According to the experiments in Section 6.2.1 in the paper, it can be seen that there are strong conflicts in some objective functions, which generate a large number of NAN values and illegal values so that the convergence of all objective functions cannot be achieved.
Since NSE, LNNSE, R E , and R 2 are the objective functions for finding the maximum, other objective functions are optimized for the minimum values. Therefore, to unify the types of the extreme value of these objective functions and to easily compare the values of objective functions, these four objective functions have been modified and converted into minimum types which are shown in equations (37)-(40):

Selection and Evaluation Framework.
To comprehensively evaluate these objective functions and their combinations which are summarized in Table 2, we constructed a simple selection and evaluation framework of objective functions which is shown in Figure 1. Firstly, a library of objective function combinations for hydrological model parameter optimization is established. Assuming that there are K kinds of combinations of objective functions, the library can be shown in the following equation: where OFC i is the i-th combination of objective functions, i � 1, 2, . . . , K, and j � 1, 2, . . . , Li is the number of the objective functions in the i-th combination. The objective functions in each combination can be partially identical but cannot be fully identical; thus, Li is different for OFC i . Then, the RMOABC optimization algorithm is employed to solve the parameter optimization problem of Xinanjiang model one by one for the objective function combination, and the corresponding Pareto optimal solutions (i.e., parameters sets) for Xinanjiang model are obtained. According to the model parameters in the Pareto solution, other remaining objective function values are calculated one by one. Finally, the relative merits of the optimization scheme corresponding to the Pareto solution obtained for each objective function combination are compared using the method of entropy-based TOPSIS.
Take parameters in order from the obtained Xinanjiang model parameter sets. Evaluate and analyze the results of the combination of objective functions.
...  Computational Intelligence and Neuroscience The selection and evaluation framework of combinations of objective functions is a simple method for comparing the optimization performance of different objectives combination, and the specific steps are described as follows: (1) Select one combination in turn from the combinations library of objective functions which is defined as (2) Use the RMOABC optimization algorithm to call the Xinanjiang model for simulating the basin runoff of the Heihe River. The noninferior solution sets (optimal model parameters) under all combinations of objective functions can be obtained, as shown in the following equation: where i = 1, 2, . . ., m; j = 1, 2, . . ., ni; and OtherSolution_OFC i are other objective functions to the i-th combination of the objective function.
where i = 1, 2, . . ., m; j = 1, 2, . . ., ni; and AllSolution OFC i is the noninferior solution set obtained by the i-th combination of objective functions. (6) Establish a decision matrix based on the entropybased TOPSIS method to sort all noninferior solutions of the combination of objective functions. (7) Comprehensively evaluate the objective function combination. According to the TOPSIS ranking results of the noninferior solution sets obtained in the foregoing operations, quantitatively describe and compare the influence of each combination of objective functions on the parameters optimization of the hydrological model.

Experiment Preparations.
The experiments of the simulation were taken on the multiobjective evolutionary algorithm (MOEA) framework [57] in the paper. The MOEA is a powerful and efficient platform and provides a Java library of open source for developing and evaluating the multiobjective optimization algorithms. The hardware environment of simulation experiments is a PC with 4 cores of Intel i7 2.6 GHz CPU and 8 GB of RAM.

Parameter Settings.
Because the selection of algorithm parameters can greatly affect the execution performance, we adopted the recommended parameter settings from the related research studies [12] for the RMOABC algorithm. The swarm size of bees (NP) is set to 100, the external archive capacity is set to 100, and the adaptive grid number is set to 25. For the Limit parameter, it is set to 0.25 * NP * D [58]. D denotes the dimension of decision variables of the optimized problem, that is, the number of Xinanjiang model parameters, 10. The dimension of the objective functions is set to the number of objective functions in the selected combination. Based on our previous experience [12], for the calibration of the parameters of Xinanjiang hydrological model, a maximum number of model evaluations in the range 2000-4000 normally ensure an efficient calibration. Thus, the stopping criterion for the optimization algorithm is the maximum number of model evaluations, and it is set to 3000.

Evaluation Indicators of Multiobjective Optimization
Algorithms. An important issue in multiobjective optimization is to evaluate and compare the quality of solution sets. The straightforward way to compare the quality of solution sets is visualization, but visual comparison cannot quantify the difference between solution sets and also becomes harder with more objectives involved [59]. Thus, quantitative evaluation methods can be interpreted as how well it represents the Pareto front and used to measure the performances of the multiobjective optimization algorithms [28]. In general, the quality indicators of solution set can be divided into four aspects: convergence, spread, uniformity, and cardinality. Convergence of a solution set refers to the closeness of the set to the Pareto front which denotes the Computational Intelligence and Neuroscience distance between the calculated noninferior front and the known true noninferior front or approximate true noninferior front. Spread of a solution set considers the region of the set covering. Uniformity of a set refers to how even the solution distribution is in the set, and an equidistant spacing amongst solutions is desirable. Cardinality of a solution set refers to the number of solutions in the set [59].
According to the literature [59], the most frequently used evaluation indicators were selected to evaluate all the four aspects of performance of the nondominated solution sets in this paper. The generational distance (GD) [60] was selected to evaluate the convergence performance, and the spacing (SP) [61] was selected to evaluate the uniformity performance. The inverted generational distance (IGD) [62] was chosen to evaluate the comprehensive quality of convergence and spread, and the hypervolume (HV) [63] was also chosen to evaluate the comprehensive quality of convergence, spread, and cardinality. Due to the desirable practical usability and theoretical properties, the HV indicator is arguably the most commonly used quality indicator and is suitable for many real-world optimization scenarios. And the HV indicator is the only popular unary indicator having the strict Pareto compliance which implies that only the Pareto front achieves a unique optimal value for a problem [59]. We also added the computational time (Time) which was recorded in seconds to evaluate the execution efficiency of the algorithms.

Model Parameter Optimization
Result. According to the above experimental scheme and algorithm parameter setting, the RMOABC optimization algorithm is used to optimize the parameters of the Xinanjiang model based on the nine objective function combinations of the hydrological model parameters compiled in Section 5.2. Each group of experiments was carried out 30 times, and the optimal parameter set of each group of experiments was collected together, and the dominant solution in each group was removed. As shown in Table 3, a total of 8,435 optimal parameter solution sets are obtained. The Pareto optimal solutions obtained by the objective function combination 1, 3, 4, 5, and 8 are relatively small. They are 95, 120, 178, 325, and 355, respectively. The objective function combination 6, 7, and 9 obtained more than 1,000 Pareto optimal solutions, and in particular, the objective function combination 2 obtained a total of 3,629.
However, since the nine objective function combinations include 19 objective functions, we take each parameter solution from each objective function combination in turn according to the steps of the selection and evaluation framework and then calculate the corresponding values of the remaining objective functions. For example, in the objective function combination 2, four objective functions are included, and one parameter solution is sequentially extracted from the corresponding 3,629 solutions, and the runoff data are simulated by the Xinanjiang model, and then the remaining 15 objective functions' values are calculated.
In this way, we can get all the values of 19 objective functions corresponding to the 8,435 optimal parameter solution sets.
Nevertheless, we found that there are many NAN values in the newly calculated objective function values which are mainly concentrated in the objective functions of LNNSE, E MSL , E TRMS , E MSR , and E MA L . After analysis, the main reasons for these NAN values are as follows. (1) As there is a value of 0 in the observation data or the prediction data, the natural logarithm calculation in LNNSE and E MSL will generate the NAN value. (2) In the hydrological model, 0 is taken as the divisor. The occurrence of the NAN value in the objective function affects the parameter solution which cannot be sorted. Therefore, we removed all the solutions which contain NAN value in the objective function. After elimination of all the solutions containing NAN values, a total of 6,918 remaining multiobjective parameter solutions are obtained.
As described in Section 5.2, in order to facilitate the comparison of these 19 objective functions, we have converted all of these objective functions into minimization functions, and the function values should be in the range [0, 1], and the closer the function value is to 0, the higher the degree of optimization of the solution. However, after reviewing the result set of the objective function values, it was found that a large amount of data larger than 1 exists in the remaining objective function values. It shows that there are strong conflicts between these objective functions. The optimal parameter set obtained for some objective functions does not perform well or even worse on other objective functions. To compare the performance of the objective function, we removed all model parameter solutions containing objective function values greater than 1 from the experimental result set. After eliminating the solutions with NAN value and value greater than 1, a total of 165 optimal parameter solution sets are obtained, and there are about 6.900 solutions whose objective function value exceeds 1.
Next, we compared the Pareto dominance of 19 objective function values with 165 optimal parameter solutions and found that there is no dominant solution in these optimal parameter solutions. Therefore, the nine objective function combinations finally obtained 165 optimal parameter solutions.
The number of optimal solution sets of each objective function combination in the total number of solution sets is shown in Table 4 and Figure 2. It can be seen that the objective function combination 2 and the objective function 7 obtain more optimal solutions, respectively. Their proportions are 45.45% and 26.67%, respectively. The optimal solutions obtained by the remaining objective combinations are relatively small, especially the objective function combinations 3 and 4 which only obtained two optimal solution sets. The objective function combinations 2 and 7 also obtained 3,629 and 1,667 solution sets before the data preprocessing, respectively, which indicates that these two objective function combinations have stronger ability to obtain the optimal parameter solution sets, and the ability of their obtained parameter solution sets is also strong for adapting to other objectives.

Calculating the Closeness of the Optimal Solution Sets of Objective Function Combinations.
After the above preprocessing for the Pareto optimal parameter solutions generated by the RMOABC algorithms based on the 19 objective functions for a validation period, a total of 165 model optimal parameter sets are obtained for the Xinanjiang model. The sample Pareto optimal parameter solutions are shown in Table 5.
The Pareto optimal solution of the hydrological model contains 19 optimization objectives, that is, there are 19 attribute values. Except the optimal solution corresponding to the attribute value of the objective function    Then, the decision matrix D will be solved according to the steps of the entropy-based TOPSIS method described in Section 5.1. In the TOPSIS method, the weight of each objective function is the key factor affecting the ranking result of optimal solutions, but the common method is to assign the corresponding weight values to the objective functions based on the experience of the decision maker for the specific problem. In this paper, 19 objective functions have been selected for the research; therefore, we cannot assign weights to these objective functions reasonably and scientifically. Thus, we adopt the method of calculating the information entropy of each objective function in the optimal solutions to determine their weights for more accuracy [38]. Table 6 indicates the entropy value and weights of 19 functions which were calculated according to equations (29) and (31). From the table, it can be found that the 19 objective functions show unequal weights. R E function shows the highest weight, and NSE and LNNSE also have an important effect on the model parameters optimization. E ROC indicates the lowest weight, 0.00000, which means the objective function almost has no effect on the hydrological model optimization results. And E M4 , E MSD , E SFDC , and E MSR also have little importance for the optimization problem.
After the completion of the weight calculation for each objection function, the weighted normalized decision matrix Z can be calculated. Then, the positive ideal solution (PIS) and negative ideal solution (NIS) will be determined. And the separation measures by the m-dimensional Euclidean distance for each alternative will be obtained. Finally, based on the entropy-TOPSIS method, the overall performance coefficient of each nondominated solution which is expressed as relative closeness to the ideal solution will be obtained according to equation (36). Since the number of optimal solutions is large, Table 7 shows all the 19 objective function values of the sample Pareto optimal parameter solutions in Table 5 and their relative closeness which is calculated by the entropy-based TOPSIS method.

Ranking the TOPSIS Preference Order.
Finally, by sorting the closeness coefficients C * i in descending order, the ranking order of all alternatives (i.e., Pareto optimal solutions) can be obtained. The higher value of C * i means the higher priority; the smaller value of C * i means the lower priority. The maximum value of the overall ranking of closeness coefficients calculated by the 165 optimal parameter solution sets is 0.95196 which means the solution is the optimal alternative, and the minimum value is 0.25942 which means the related solution is the worse alternative. To compare the optimal solution sets obtained by each objective function combination, we multiply the value of closeness coefficients by 100 and convert it into a 5-point system, that is, the converted value is "Excellent" between 90 and 100, "Good" between the 80 to 89, "Medium" between 70 and 79, "Pass" between 60 and 69, and "Poor" for values lower than 60.
After dividing the ranking score into 5 levels, we then analyzed the distribution of the optimal parameter solution sets obtained by the nine combinations of the objective functions in the 5-segment position interval. The results are shown in Table 8 and Figure 3. We can see that in the 165 optimal solution sets, objective combinations 2 and 7 obtain a larger number of optimal solutions, while other objective combinations obtain fewer optimal solutions. Taking objective combination 1 as an example, five of its six optimal solutions are above the "Pass" level, indicating that the quality of the solution corresponding to the objective combination is higher, but its number of solutions is very small, accounting for only 3.64% of the total solution. It shows that the ability of these objective function combinations to obtain the optimal parameter solution sets is not strong, and the obtained parameter solution set cannot adapt to other objective function requirements which means that they do not have advantages in the competition.
It can be seen from the distribution diagram of TOPSIS sorting position in Figure 3 that the number of distributions of objective combinations 2 and 7 in the five sorting intervals is very uniform, especially in the Excellent and Good intervals, and the solutions of the top ranking are superior in number and ratio to the other objective functions. It shows that objective function combinations 2 and 7 can achieve more and better solutions than other objective combinations.
To compare the quality of the solutions corresponding to each objective combination, we extracted all the solutions that reach the excellent level and sorted them according to the closeness value. A total of 17 Pareto optimal parameter sets have reached the excellent level whose relative closeness values are bigger than 0.90, and the ranking results are shown in Table 9. The solution having the highest TOPSIS relative closeness indicated with an asterisk in the table is the compromise solution which fits the decision maker's preference most.
It can be seen from the table that among the top 10 solutions, the first to fourth and the sixth are the optimal solution sets corresponding to the objective combination 2,

16
Computational Intelligence and Neuroscience and in all 17 solution sets, the objective combination 2 has 6 which accounts for 35.29%. It shows that the optimal solution corresponding to the objective combination 2 not only has higher quality but also has stronger competitiveness and can basically adapt to the goals and requirements of all 19 objective functions. The other outstanding objective function combination, i.e., combination 7, also has 6 solution sets to reach the excellent level, but its corresponding Table 8: The distribution of the optimal parameter solution sets obtained by the nine combinations of the objective functions in the 5segment position interval.

ID Number Excellent
Good Medium  Pass  Poor  Sum  1  0  3  2  0  1  6  2  6  19  14  18  18  75  3  2  0  0  0  0  2  4  0  1  0  1  0  2  5  0  1  1  0  1  3  6  0  6  3  2  1  12  7  6  7  14  11  6  44  8  2  3  1  2  3   highest order solution is only in the 7th position, and the rest of the solutions are ranked in the lower positions. The solutions related to the remaining objective combinations do not perform very well. For example, the objective combinations 3 and 8 only have 2 solutions, and the objective combination 9 only has 1 solution in the excellent level. Not only the quality of the solution is low, but the number of solutions is also very small. And for the objective combinations 1, 4, 5, and 6, they even do not have the optimal solution set to achieve the "Excellent" level. Therefore, it can be seen from the experimental results that the solutions corresponding to the objective function 2 and the objective function 7 are more widely distributed than the other objective combinations, and their performances are better, especially the solution quality corresponding to the objective function 2 is more excellent.

Comparison of Optimal Parameter Results
Corresponding to Each Objective Function Combination. First, each parameter set is named in turn by the sequence number sorted by TOPSIS. Then, according to the ID number of each objective function combination, the solution set with the highest ranking of each target function combination is extracted from the TOPSIS sorted parameter solution set, and the 10 optimized parameter values, the corresponding closeness values, and the corresponding sorting results are shown in Table 10. For example, the objective combination 2 corresponds to solution 1, the objective combination 8 corresponds to solution 5, and so on. Meanwhile, the 19 objective function values for optimal parameter solution related to the 9 combinations of objective functions are listed in Table 11. The solutions of the objective combinations 2, 8, 7, 9, and 3 are at the "Excellent" level, and the solutions corresponding to the remaining objective combinations are at the "Good" level. In Section 5.2, we made corresponding modifications to the objective functions and converted the maximum objective function to the minimum type, and the values of each objective function were normalized whose value range is limited in [0, 1]. The smaller the value, the better the solution related to the objective function. Therefore, to compare the relative differences of the attributes of the nine objective function combinations corresponding to the optimal solutions, we compare the 19 objective function values in Table 11, and the comparison results are shown in Figure 4.
It can be seen from the figure that among these optimal solutions, the parameter set of the objective function combinations 2, 8, and 7 performs better, and smaller value can be obtained in each objective function. It shows that the objective combinations can meet the optimization requirements of these 19 objective functions and can achieve better optimization results. Among them, the objective function combination 2 performs particularly well, and R E and P BIAS clearly achieved better results. However, the corresponding solutions of objective combinations 9 and 3 are also in excellent grade, but they have some fluctuations in the objective function values; that is, they perform better in some objective functions but perform poorly in other objective functions. Especially for the solutions corresponding to the objective combinations 1, 6, 4, and 5 in the "Good" level, the fluctuations are more obvious. Taking the objective combination 5 as an example, the optimal optimization effect is obtained on E RMOV , but it achieved the worst optimization results in LNNSE, E MSL , R E , and P BIAS . It is indicated that the solutions related to these objective functions are not stable enough to fully cover the requirements of the 19 objective functions and they do not have strong adaptability.

Result Comparison of Multiobjective Evaluation
Indicators. It can be seen from the above experimental results that compared with other objective function combinations, the objective function combinations 2 and 7 can obtain a widely distributed and high-quality parameter solution set. In order to compare the quality of the corresponding solutions of these two objective function combinations, the generational distance (GD), spacing (SP), inverted generational distance (IGD), and Hypervolume (HV) were adopted to evaluate the performance of convergence, uniformity, spread, and cardinality of the solution set. And the computational time (Time) was used to evaluate the execution efficiency of the algorithms.
Firstly, according to the approximate noninferior solution set obtained in the previous experiments and the noninferior solution sets obtained by the 30 independent runs of algorithm, the five quality indicators of noninferior solutions under each objective combination are calculated, and then the performances of the two objective combinations on the optimization of hydrological model parameters can be quantitatively described and compared. The five performance metrics for parameter optimization problem are shown in Table 12. Mean and SD represents the average  Table 12, the conclusions consistent with the above can be obtained. (1) The optimal noninferior parameter solution sets corresponding to objective combinations 2 and 7 have good performance in convergence, uniformity, spread, and cardinality. (2) In comparison, the optimal parameter solutions of objective function combination of 2 have relatively better results for generational distance (GD), spacing (SP), inverted generational distance (IGD), and hypervolume (HV). It demonstrated that the combination 2 is better than the combination 7 in terms of the performance indicators. (3) Because the number of objective functions in combination 2 is bigger than combination 7, it has to take more execution time to obtain the final solution set. (4) Overall, the objective combinations 2 and 7 can effectively coordinate between the 19 various conflicting objectives which are arranged in the paper, and it is expected to provide forecasters with a more comprehensive basis for decision-making.

Conclusions
In the multiobjective parameters' optimization process of the hydrological model, a combination of efficiency criteria (objective functions) is defined as the mathematical measures to evaluate how well the model simulation fits the available observed data. There are a large number of objective functions, and each objective function may place more or less emphasis on different aspects of the hydrological model. However, how to choose and evaluate the different combination of objective functions efficiently has been a challenge for even the most experienced hydrologists. To help the researchers select and evaluate the objective functions that are suitable for their research, nine objective function combinations consisting of 19 objective functions were collected from the literature and arranged in this paper.
Then, a selection and evaluation framework of objective function combinations was presented. The daily runoff forecasting model, Xinanjiang model, is taken as an application case, and the RMOABC algorithm is introduced to find the nondominated parameter solutions for handling the optimization problem of the hydrological model in the Heihe River basin. And a simple information entropy-based TOPSIS ranking method is adapted to rank the nondominated parameter solutions and compare the performance of optimization results of different objective function combinations.
From the optimization results, it can be found that there are a large number of conflicts between the objective functions of objective combinations, which leads to the result that optimal parameter solutions are difficult to satisfy the optimization tasks of all objective functions. Therefore, the combination of objective functions should be chosen carefully and reasonably; otherwise, the obtained optimal parameters will not reflect the real characteristics of the hydrological model and thus cannot complete the hydrological model runoff forecasting task well.    Then, the information entropy-based TOPSIS ranking method was taken to calculate the closeness of the optimal solution sets of the objective function combinations, and the ranking order of Pareto optimal solutions was obtained by sorting the closeness coefficients C * i in descending order. Based on the analysis of the optimization results, the number of distributions of the objective combinations 2 and 7 in the five sorting intervals is very uniform, especially in the "Excellent" and "Good" levels. And their Pareto solutions with the highest ranking have a great advantage in the number and proportion of the solution sets to other objective function combinations which shows that objective function combinations 2 and 7 can achieve more and better solutions than other combinations. Then, we compared the optimal parameter results corresponding to each objective function combined with the multiobjective evaluation indicator results. The solution corresponding to objective combination 2 has better performance in the quality indicators.
The experiment results showed that the convergence and distribution of combination 2 of objective functions are better than the other combinations, and it can effectively balance the performance of various objectives. It also showed that the entropy-based TOPSIS can provide an effective and better method to calculate objective weights that permit trade-offs between attributes and rank the nondominated solutions impersonally for the parameters' optimization problem of the hydrological model, and it can also provide a more comprehensive decision basis for hydrological forecasting. The future work is to construct a more accurate and balanced combination of objective functions that has good distribution and convergence in different objective spaces.

Data Availability
The observed data for the Xinanjiang model can be downloaded from the website of Scienific Data Center of Cold and Arid Regions: http://westdc.westgis.ac.cn.