Research Article  Open Access
Parameter Estimation in RainfallRunoff Modelling Using Distributed Versions of Particle Swarm Optimization Algorithm
Abstract
The presented paper provides the analysis of selected versions of the particle swarm optimization (PSO) algorithm. The tested versions of the PSO were combined with the shuffling mechanism, which splits the model population into complexes and performs distributed PSO optimization. One of them is a new proposed PSO modification, APartW, which enhances the global exploration and local exploitation in the parametric space during the optimization process through the new updating mechanism applied on the PSO inertia weight. The performances of four selected PSO methods were tested on 11 benchmark optimization problems, which were prepared for the special session on singleobjective realparameter optimization CEC 2005. The results confirm that the tested new APartW PSO variant is comparable with other existing distributed PSO versions, AdaptW and LinTimeVarW. The distributed PSO versions were developed for finding the solution of inverse problems related to the estimation of parameters of hydrological model Bilan. The results of the case study, made on the selected set of 30 catchments obtained from MOPEX database, show that tested distributed PSO versions provide suitable estimates of Bilan model parameters and thus can be used for solving related inverse problems during the calibration process of studied water balance hydrological model.
1. Introduction
Particle swarm optimization (PSO) was established in 1995 by Kennedy and Eberhart [1]. It is an evolutionary optimization technique inspired by a behaviour of population of individuals, which form groups or swarms like flock of birds or school of fishes. It works with population of particles, which are finding the optimal solution during their search within the search space by collaboration between individuals and by exchanging information about their best position in the space.
The PSO belongs to the stochastic, metaheuristic evolutionary computational techniques, which are based on the swarm intelligence [2]. The main benefits of this optimization are the small number of parameters, which need to be adjusted, and an easy implementation. When comparing with standard local search methods, the PSO does not require a knowledge about the gradient of the optimized function [3, 4].
Its optimization process often suffers from premature convergence or from trapping in the local optimum. Therefore, the recent PSO research focuses on the development of new adaptation strategies. The most important adaptations improve the particle’s velocity estimates by the adaptation of PSO inertia weight [5, 6].
Special attention is also put on the development of distributed version of PSO [7, 8]. The original population is divided into subswarms or complexes. The complexes are simultaneously evolving over the parametric space, and they periodically exchange information according to some prescribed migration rules [9, 10].
The PSO has been successfully applied into many reallife optimization problems in engineering. In recent years, the PSO optimization significantly enhances the estimation of parameters of hydrological models [11–14].
For example, Jiang et al. [15] applied PSO for calibration of the rainfallrunoff model HIMS. They compared classical PSO algorithm with distributed PSO versions, which are using the complexes and shuffling mechanism. They found out that the distributed PSO variants are significantly better than the original one.
ZambranoBigiarini and Rojas [16] developed the standalone global optimization for calibrating the hydrological models based on the PSO (called hydroPSO and is available for R interpreter). The methods enable testing different PSO versions on calibration of analysed hydrological model. The group of tested PSO versions showed good optimization performance, and it was an effective and efficient tool for calibration of both surface and groundwater hydrological models.
The main aim of the presented paper is to test the selected PSO distributed versions on singleobjective benchmark optimization problems and to apply them on calibration of hydrological model Bilan. The case study is conducted on 30 US catchments, for which the data of hydrological and meteorological forcings are obtained from MOPEX experiment [17].
The rest of the paper is organized as follows. Section 2 provides details of the standard PSO algorithm and its tested modifications. Section 3 explains methodology used during singleobjective benchmark optimization problems and methodology of the rainfallrunoff model simulations with description of the Bilan rainfallrunoff model. Results are summarized in Section 4 and discussion is provided in Section 5. Finally, the main findings are concluded in Section 6.
2. Particle Swarm Optimization
This section provides description of tested versions of distributed PSO. At first, the original PSO is described, and then the analysed variants of PSO are provided. The distribution strategy is explained in the last section.
2.1. The Original PSO
The original PSO algorithm estimates a new particle’s location using information of particle’s velocity. The velocity is updated as and the location is simply defined as for all , where is the total number of particles in swarm population, , where is the total number of dimensions (i.e., the number of parameters of hydrological model), and denote acceleration constants predefined by the user, and and are independent random vectors sampled from a uniform distribution in the range .
The component with particle’s previous best position in (1) represents the cognition knowledge of swarm particles. When optimization problem is the minimization, for all so far known locations of given particle . The component with the best position among all particles in (1) controls the social influence of swarm particles. For all in the population . is the analysed singleobjective function [1, 19].
The particle’s locations are in our research initialized randomly within the search space using the initialization based on the Latin hypercube sampling technique [20, 21]. Particle’s velocity is initialized randomly from the interval , where is equal to [22, 23].
2.2. Analysed Modifications of PSO
We analysed four versions of the PSO algorithm. They differ according to the applied particle’s velocity adaptation. All versions were tested as asynchronous distributed PSO. The modifications are the following.
ConstrFactor. In the first modification, the parameter of constriction factor is implemented into the PSO algorithm [24]. If the adapted particle’s velocity with is as then where and .
LinTimeVarW. In the second PSO modification, the linearly decreasing inertia weight is used [25]. The adapted particle’s velocity equation is shown in The value of inertia weight decreases at each iteration linearly from to and is simply defined as
AdaptW. One of the adaptive inertia weights from [26] was used in this paper. The modification of the inertia weight parameter is as for which where shows the success of th particle, is the success percentage of the swarm, is the size of the population, , and . The adaptive inertia weight value is updated based on a feedback parameter which is in this case the variable of percentage of success.
The AdaptW PSO method provided the best performance in the original paper of Nickabadi et al. [26], and it was also the best modification of inertia weight out of total distributed variants of PSO (for details see [27]).
APartW. The proposed new adaptive inertia weight modification is based on the current position of the particle, and it combines the global exploration and local exploitation in the space. The value of the inertia weight parameter is updated at each generation according to the development of particle’s location. If the particle improves its position compared to the best location achieved so far, it is closer to the searched optimal value. So, if , the local exploitation is supported. The inertia weight is calculated as and thus, for random vector sampled from a uniform distribution in the range , and ; the resulted inertia weight is between .
On the other hand, if the particle does not achieve better position, the global exploration is encouraged. So, if , the resulted inertia weight is computed as where the variables are the same as in (9) and the inertia weight is thus between .
2.3. Distributed Version of PSO
During the optimization process, a premature convergence to a local optimum could appear. Due to this, different distribution strategies of the swarm were developed [27, 28]. A new distribution strategy of the swarm called shuffled complex evolution (SCE) was proposed by [9]. In this paper, the SCEPSO method introduced by [29] is used, where the shuffled complex evolution is combined with the PSO optimization. The only difference in our research is that all particles from the complex are participating in the PSO algorithm, not only a predefined number of them.
In the SCE strategy, after initialization of particle’s position and velocity, the whole population is divided into a predefined number of complexes . The division is according to the functional value of each particle. All particles are sorted in increasing order, and then each th complex receives particles [30]. Each complex simultaneously searches through the parametric space, and, after a predefined number of iterations in one complex, all particles return to the swarm. The shuffling of particles and redistribution into the complexes are made, and the process is repeated until the termination criteria are satisfied.
The shuffling mechanism preserves the population diversity and helps to prevent premature convergence. The SCEPSO was already applied for comparison of modifications of PSO in our previous research [27], and it was found that the SCE strategy gives significantly better results.
In this paper, all four PSO variants were extended into a distributed version using SCEPSO technique. Since the APartW algorithm is a new proposed version, the simplified PSO algorithm using APartW adaptation and SCE mechanism is shown in Algorithm 1, where is the number of particles in one complex. The other PSO modifications differ only on lines 9, 11, and 16 in Algorithm 1.

3. Methodology of SingleObjective Optimization Problems
We analysed the tested distributed versions of PSO on two sets of singleobjective optimization problems. All optimization problems were minimizations. The total number of analysed optimization problems was .
The first set is represented by the benchmark problems, which were specially prepared for CEC 2005 singleobjective optimization session [31]. The total number of optimization runs was (i.e., benchmark functions × 4 PSO variants × 25 program runs).
The second set consists of optimization problems. On datasets of MOPEX catchments, we analysed benchmark questions, which are standard objective functions used for solving inverse problem related to calibrations of hydrological models [32, 33]. Each optimization based on one objective function and one data forcing for one catchment was repeated 25 times. The total number of optimization runs was 12 000 (i.e., 4 objective functions × 30 catchments × 4 PSO variants × 25 program runs).
All algorithms, benchmark functions, and Bilan model were written in C++ programming language, and the code ran under 64bit Linux operating system. All graphs and statistical tests were made in R statistical software environment [34].
3.1. The CEC 2005 Benchmark Optimization Problems
In order to compare the optimization algorithms, the singleobjective benchmark optimization problems are analysed. We used 11 benchmark functions from the special session on single optimization problems CEC 2005 [31]. This set of benchmark functions was used by many researches for solving optimization problems [27, 35, 36].
We tested the following minimization benchmark problems: sphere (), Schwefel 1.2 (), elliptic rotated (), Schwefel 1.2 with noise (), Schwefel 2.6 (), Rosenbrock (), Griewank rotated (), Ackley rotated (), Rastrigin (), Rastrigin rotated (), and Weierstrass rotated (). The optimization problems are unimodal functions; the remaining are multimodal functions. For formulas of all functions, range of the search space, location of the minimum value, and other descriptions, see [31].
The setting of distributed PSO versions follows [22, 27, 37]. The number of complexes is set to 6, and there are 25 particles at each complex. The number of shuffling is , and the maximum number of function evaluation is set to . The problem space has dimensions due to preservation the setup from [27]. For results analysis, the total number of optimization runs is set to , where each run stops when the maximum number of function evaluations is achieved. In terms of the PSO coefficients, the acceleration constants for modifications with inertia weight are , and the acceleration constants for adaptation using constriction factor are with constriction coefficient .
3.2. Optimization of the Hydrological Model: Case Study
After comparison of the analysed optimization algorithms on benchmark functions, the developed distributed versions of PSO were applied on solving the reallife optimization related to the estimation of values of parameters of lumped hydrological model Bilan.
The settings of PSO parameters on Bilan optimization problems were selected after running several tests with different parameter settings. The number of complexes is set to . Each complex is composed of 40 particles, and the number of generations in one complex is 20. The number of shuffling is . To analyse the results, the total number of model runs is . The acceleration constants and value for constriction factor are the same as in the singleobjective benchmark optimization problems.
3.2.1. The Hydrological Model Bilan
In the case study, the calibration of Bilan rainfallrunoff model is studied. Bilan is a lumped physically based water balance model developed in T. G. Masaryk Water Research Institute in the Czech Republic [38]. It is a standard tool commonly used for assessment of water balance in the catchment [39–41].
It is a conceptual hydrological model, which explains the hydrological balance of a catchment using the system of mathematical relationships, which preserve mass balance. It describes basic principles of water balance on ground, in the zone of aeration, including the effect of vegetation cover and groundwater. The main model forcings are time series of precipitation [mm], air temperature [], and relative air humidity . Time series of air temperature and relative air humidity are used to estimate the potential evapotranspiration [39, 42].
The temporal dynamics of reservoirs is described by firstorder differential equations, which are numerically solved using the Euler method. The calculated total streamflow is given by two components of the river flow. The fast response is simulated through direct runoff reservoir, and the second slow runoff component is explained with baseflow reservoir [43]. The Bilan scheme is shown in Figure 1, and further description of the model could be found in [42].
The total amount of parameters of the daily version of the Bilan model is six, and they are listed in Table 1. The search space was constrained with physically meaningful ranges of parameters (see column PC in Table 1). The parameter constraints were estimated by an expert knowledge.
 
PC = parameter constraints. 
3.2.2. Objective Functions
The Bilan model is calibrated against the observed streamflow data, so the model time series of observed streamflow are used for calibrating of the model. Therefore, different calibration indices based on information obtained from model residuals are used for estimation of Bilan parameters. The solution of related inverse problem, which minimizes the analysed hydrological index, objective function, was used. This approach is a standard way of calibration of lumped hydrological models [44, 45].
The investigated objective functions are in hydrological modelling commonly used accuracy criteria: mean squared error (MSE), mean absolute error (MAE), mean absolute percentage error (MAPE), and NashSutcliffe efficiency (NS) [32, 33, 46].
The analysed objective functions are defined as where is observed streamflow [mm], is modelled streamflow [mm], is mean of the observed streamflow [mm], and is the total number of observations.
The fraction residual variance () was used for the calibration purposes; however, its extension the NashSutcliffe coefficient is used for standardized assessment of the results of calibration of hydrological models Nash [47, 48]. is defined as
3.2.3. Dataset
For evaluation of the optimization ability of the proposed PSO modifications, we calibrated Bilan using datasets from US catchments. The meteorological and hydrological data were obtained from Model Parameter Estimation Experiment project (MOPEX) [17], which serves for benchmarking of hydrological models and calibration approaches.
For the analysis, the daily records from to were used. The main meteorological forcings of Bilan model were mean areal precipitation, mean air temperature calculated as an average of given maximum and minimum daily temperature, and potential evaporation. Table 2 lists the major characteristics of each catchment including latitude, longitude, drainage area, dominant soil, and vegetation types.
 
MF = mixed forest; CS = closed shrublands. 
4. Results
4.1. The CEC 2005 Benchmark Optimization Problems
Figure 2 presents the convergence graphs for each benchmark function while utilizing all four distributed modifications of the PSO algorithm. The axis indicates the number of function evaluations and the axis the logarithmic value of the best fitness, that is, the difference between the searched and the best achieved functional value. A decline with the number of evaluations is clearly visible for LinTimeVarW, AdaptW, and APartW modifications in all functions, and it thus indicates the approach to the global optimum. The ConstrFactor variant has the worst performance, where the decline towards the global optimum is not evident.
For the statistical comparison, the nonparametric Wilcoxon test was used according to [49]. Inputs to the testing procedure were the best fitness values achieved for the PSO modifications. The optimization performance of the new proposed APartW was statistically compared only with the AdaptW variant because the AdaptW was the best modification in the research papers of [26, 27], and thus it serves as a benchmark method.
Table 3 summarizes results of the Wilcoxon test, and other statistical indices are also indicated. The minimum, 25% quartile, median, 75% quartile, maximum, mean, standard deviation, and value are available. All program runs of each modification were compared in each numbered evaluation. The values reflect the best achieved fitness and report other statistical indices belonging to the same numbered evaluation.
 
value of the statistical Wilcoxon test. Nonsignificant differences in medians are highlighted in bold; significantly better median in terms of the new adaptive PSO APartW is highlighted in bold with stars ( for value ≤ 0.05; for value ≤ 0.001). 
The results of the analysis in Table 3 show that the APartW modification gives significantly better results for three benchmark functions (, , and ). In functions , , and , there is no significant difference between APartW and AdaptW. Beyond that, in functions and , both APartW and AdaptW found the global minimum. For multimodal functions , , and , the AdaptW variant gives significantly better results. The differences in the APartW and AdaptW are in agreement with the convergence graphs on Figure 2.
Based on our findings, we can conclude that the APartW version of PSO is comparable with AdaptW modification in optimizing chosen benchmark functions. The APartW is reliable for singleobjective optimization and thus suitable for calibration of the Bilan rainfallrunoff model parameters, where only one objective function at time is minimized.
4.2. The Calibration of Bilan Model
For calibration of the Bilan model, hydrological and meteorological data from 30 US catchments were used. The analysis of the main findings from all catchments is provided, but for clarity, the detailed results from only one catchment are displayed. The chosen catchment (USGS ID ) provides the outline of obtained results, and it serves to perform the overall findings.
We used the ensemble simulations of the total number of model runs equal to . At each ensemble simulation, only one objective function was optimized. Within all simulations, we also calculated the remaining three accuracy criteria for evaluating the overall performance.
The results in Table 4 show the best achieved values of objective functions within all catchments. The optimized objective function is in the first column; the calculated criteria are in the last four columns. It is evident that the APartW method provides similar results to the LinTimeVarW version. They both reached the best values in four cases (highlighted in bold in the table). In addition, the ConstrFactor gave the worst results (reached the best values in two cases), and, on the other hand, the AdaptW method achieved the best results (reached the best values in six cases). Even though the APartW method does not provide the best value of the optimized objective function, it achieves good results for the other criteria which are not optimized.
