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.
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.
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.
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.
In order to determine whether the implementation of the distributed PSO modifications into Bilan model increases the model performance, we compared our results with the integrated binary search (BS) optimization technique. The best achieved objective criteria for all catchments using BS method are the following: , , , and . The results are for optimization based on , , , and , respectively. In comparison with Table 4, the optimization with PSO method is always better than with the integrated BS method.
Table 5 displays results from the contrast test of the unadjusted medians according to [49]. After pairwise comparison of all PSO modifications, we determined the ranks of each method. The APartW variant achieved the first rank once, the second rank two times, and the third rank once. The best method seems to be the AdaptW, which achieved the best results two times and the second rank also two times. On the other hand, the worst is the ConstrFactor version, which was always worse than the others. Additionally, differences in medians between LinTimeVarW, AdaptW, and APartW are very small, which indicates similar performances.
In addition to the contrast test, we also conducted the Wilcoxon pair test of medians according to [49]. The ranks are displayed in the last column in Table 5. The obtained results confirm the results from the contrast test. The differences in the ranks are in the simulations based on and objective functions, where APartW variant is better than the AdaptW or as good as AdaptW, respectively. In terms of Wilcoxon test, the APartW is the best modification and ConstrFactor is again the worst.
Since the APartW modification is a new proposed PSO variant, in Figure 3 is displayed the time series of observed and modelled runoff using this method. For clarity, only one chosen year is depicted. The figure gives an example of ensemble simulations with the Bilan model, where the results from the total model runs are coloured in grey. It is evident that the envelope curve of the ensemble simulations would cover most of the observed data points. In the figure, also the streamflow calculated by the best model is plotted (red line), that is, the simulation with the highest value of .
From Figure 3 it can be seen that the simulation by the best model underestimates the runoff. This behaviour is also clear from the scatter plot displayed in Figure 4, where the regression line lies below the line. The corresponding residuals, that is, differences between observed and simulated streamflow, are depicted in Figure 5. The range of the residuals is approximately for ensemble runs and for the best model.
The values of parameters obtained by the best model of each PSO modification for the chosen catchment are in Table 6. The obtained NashSutcliffe efficiency values are also indicated in the table. It is seen that the are very similar for all PSO modifications. The APartW variant achieved the smallest value, but it was not significantly smaller than the others.
The results showed that the PSO versions using adaptive inertia weight for calibration of the Bilan model give better results than other PSO modifications. This is in agreement with Nickabadi et al. [26]. We found out that the linearly decreasing inertia weight for updating the particle’s velocity is more promising than the constriction factor, which is in contradiction with Eberhart and Shi [22], but, in our research, the distributed versions were used.
5. Discussion
5.1. The CEC 2005 Benchmark Optimization Problems
Both AdaptW and APartW modifications obtained the global minimum, where the error value was less than , in three functions. These functions are , (unimodal), and (multimodal), and the results are comparable with our previous research [27].
Better results were obtained by Liang and Suganthan [35], who applied distributed PSO algorithm with local search. They achieved the global optimum in functions , , , , and . The swarm was dynamic with frequent regrouping and the swarm’s size was very small.
Bui et al. [36] achieved the global optimum in the unimodal functions , , and using different versions of PSO. Their algorithms did not converge to the global minimum in any multimodal function, and, thus, we can consider our distributed PSO versions better.
Our results are similar to the results obtained by Nickabadi et al. [26]. Their algorithms achieved the global optimum in functions , , and .
When comparing our results with other optimization techniques within the CEC 2005 benchmark problems, our distributed PSO versions give comparable results as CoEVO algorithm [50]. They also obtained the global optimum in functions , , and . Our results are better than the DE algorithm of Rönkkönen et al. [51], which achieved the global optimum in functions and or better than BLX_MA method of Molina et al. [52], which converged only in functions and .
5.2. The Calibration of Bilan Model
Our best obtained values (>0.74, see Table 4) are comparable to a work of Brauer et al. [53, 54], who implemented the hydroPSO package [16] for calibration of the Wageningen lowland runoff simulator (WALRUS). The WALRUS is a rainfallrunoff model often used in sloping lowland catchments. The differences in their work are that, during the calibration of the model, they used year of discharge observations and optimized parameters. Their best achieved values for two catchments were and . In our work, we used years of observation and calibrated parameters for catchments.
Similar values of criteria were also obtained by Jiang et al. [55] in the calibration of the Xin’anjiang hydrological model, by Lawrence et al. [56] in calibration of the HBV model or by Zhang et al. [57] in calibration of SWAT hydrological model. Thus, we can consider our optimized parameters of the Bilan model as sufficient estimation, and therefore the used PSO distributed versions are suitable for the calibration.
The Bilan model was extended by a shuffled complex differential evolution (SCDE) global optimization method and the model was applied on catchments in the Czech Republic. The best achieved values during calibration of the Bilan model using SCDE were between and , whereas using existing integrated local optimization technique with expert knowledge it was [58]. Our obtained results were comparable with this research; however, different dataset was used.
To determine if our distributed versions of PSO improve the model performance, we compared our results with the binary search method. The binary search is a default optimization technique integrated in the Bilan model. We found out that the objective criteria obtained from the PSO optimization gave better results. Therefore, the implementation of the PSO into the Bilan model improves the model simulations.
6. Conclusions
The main aim of this paper was to test selected PSO distributed versions on singleobjective benchmark optimization problems and to apply them on calibration of hydrological model Bilan. For all PSO versions, optimization problems were analysed, in which minimizations for benchmark problems (i.e., benchmark functions × 25 program runs) and inverse hydrological problems (i.e., objective functions × 30 catchments × 25 program runs) were solved.
The new proposed variant APartW was compared with other existing PSO modifications—ConstrFactor, LinTimeVarW, and AdaptW on benchmark functions prepared for the special session on realparameter optimization of CEC 2005. The APartW version is comparable with the AdaptW and LinTimeVarW variants, whereas the ConstrFactor had the worst performance.
The statistical comparison of AdaptW and APartW modifications showed that both methods obtained the global minimum in three functions (, , and ). The APartW variant gave significantly better results in functions , , and . Beyond that, in functions , , and , there was no significant difference between the APartW and AdaptW methods.
All four PSO modifications were implemented into the Bilan rainfallrunoff model for solving inverse hydrological problems. Based on the contrast test of the unadjusted medians and Wilcoxon test, it was found out that the APartW and AdaptW variants provided the best results. The ConstrFactor performed the worst.
Our results highlighted that distributed versions of PSO are promising in singleobjective optimization problems. It was confirmed that adaptive variants of the inertia weight are better than linearly decreasing weight. It was also found out that the PSO modifications with parameters of inertia weight give significantly better results than the variant with constriction factor.
The results of this paper extended the range of utilization of the PSO global optimization techniques. The performance of the distributed versions of PSO is promising, and the algorithms can be implemented into other real optimization problems. For future work, we would like to incorporate the PSO with its distributed modifications into other hydrological models for rainfallrunoff simulations.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors are grateful to the anonymous referees for their valuable comments and suggestions for improving the presentation of this paper. This work was supported by the Internal Grant Agency, Faculty of Environmental Sciences, Czech University of Life Sciences Prague (Projects nos. 422001312315820144227 and 422001312314420154219).