Parameter Estimation in Rainfall-Runoff Modelling Using Distributed Versions of Particle Swarm Optimization Algorithm

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 shufflingmechanism, 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 single-objective real-parameter 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 fromMOPEX 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.


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][12][13][14].
For example, Jiang et al. [15] applied PSO for calibration of the rainfall-runoff 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.
Zambrano-Bigiarini 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 single-objective 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 single-objective benchmark optimization problems and methodology of the rainfall-runoff model simulations with description of the Bilan rainfall-runoff model.Results are summarized in Section 4 and discussion is provided in Section 5. Finally, the main findings are concluded in Section 6.

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.

The Original PSO.
The original PSO algorithm estimates a new particle's location using information of particle's velocity.The velocity ⃗   = (V 1 , V 2 , . . ., V Dim ) is updated as and the location ⃗   = ( 1 ,  2 , . . .,  Dim ) is simply defined as for all  = 1, . . ., , where  is the total number of particles in swarm population,  = 1, .

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  =  1 +  2 and  > 4.
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  max = 0.9 to  min = 0.4 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, Require: NC, S comp, termination criteria (1) initialize the position and velocity of all particles (2) repeat (3) E ← sorted X in increasing order according to the functional value (4) for each complex  = 1 to NC do (5) A j ← divided E into NC complexes (6) for each particle  = 1 to S comp in A j do (7) if (12) end if (13) end for (14) for each particle  = 1 to S comp in A j do (15) end for (19) end for (20) until termination criteria is met Algorithm 1: PSO algorithm with APartW adaptation and SCE mechanism.
min = 0, and  max = 1.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 27 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 ( ⃗    ) < ( ⃗  −1  ), the local exploitation is supported.The inertia weight is calculated as and thus, for random vector ⃗  sampled from a uniform distribution in the range [0, 1],  min = 0.1 and  max = 0.9; the resulted inertia weight is between [0.1, 0.5].
On the other hand, if the particle does not achieve better position, the global exploration is encouraged.So, if ( ⃗    ) ≥ ( ⃗  −1  ), the resulted inertia weight is computed as ) , (10) where the variables are the same as in (9) and the inertia weight is thus between [0.5, 0.9].

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 SCE-PSO 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 NC.The division is according to the functional value of each particle.All particles are sorted in increasing order, and then each th complex receives ⃗   , ⃗  +NC , ⃗  +2NC , . . .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 SCE-PSO was already applied for comparison of 27 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 SCE-PSO 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  comp is the number of particles in one complex.The other PSO modifications differ only on lines 9, 11, and 16 in Algorithm 1.

Methodology of Single-Objective Optimization Problems
We analysed the tested distributed versions of PSO on two sets of single-objective optimization problems.All optimization problems were minimizations.The total number of analysed optimization problems was 15.
The first set is represented by the 11 benchmark problems, which were specially prepared for CEC 2005 single-objective optimization session [31].The total number of optimization runs was 1100 (i.e., 11 benchmark functions × 4 PSO variants × 25 program runs).
The second set consists of 120 optimization problems.On 30 datasets of MOPEX catchments, we analysed 4 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 64-bit Linux operating system.All graphs and statistical tests were made in R statistical software environment [34].

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].
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 5, and the maximum number of function evaluation is set to 10 000 ⋅ Dim.The problem space has 30 dimensions due to preservation the setup from [27].For results analysis, the total number of optimization runs is set to 25, 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  1 =  2 = 2, and the acceleration constants for adaptation using constriction factor are  1 =  2 = 1.49445 with constriction coefficient  = 0.729.

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 real-life 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 3. Each complex is composed of 40 particles, and the number of generations in one complex is 20.The number of shuffling is 10.To analyse the results, the total number of model runs is 25.The acceleration constants and value for constriction factor are the same as in the single-objective benchmark optimization problems.

The Hydrological Model Bilan.
In the case study, the calibration of Bilan rainfall-runoff 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][40][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 [ ∘ C], 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.

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 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 (FRV) was used for the calibration purposes; however, its extension the Nash-Sutcliffe coefficient is used for standardized assessment of the results of calibration of hydrological models Nash [47,48].NS is defined as 3.2.3.Dataset.For evaluation of the optimization ability of the proposed PSO modifications, we calibrated Bilan using datasets from 30 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 1948 to 2003 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  the major characteristics of each catchment including latitude, longitude, drainage area, dominant soil, and vegetation types.

The CEC 2005 Benchmark Optimization Problems. Fig-
ure 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 25 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 ( 4 ,  7 , and  11 ).In functions  3 ,  5 , and  9 , there is no significant difference between APartW and AdaptW.Beyond that, in functions  1 and  2 , both APartW and AdaptW found the global minimum.For multimodal functions  6 ,  8 , and  10 , 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 single-objective optimization and thus suitable for calibration of the Bilan rainfall-runoff model parameters, where only one objective function at time is minimized.

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 01531000) 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 25.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 30 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: MSE = 1.044 + 00, MAE = 3.598 − 01, MAPE = 3.287 − 01, and NS = 7.462 − 01.The results are for optimization based on MSE, MAE, MAPE, and NS, 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 MSE and MAPE 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 25 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 NS.
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 1 : 1 line.The corresponding residuals, that is, differences between observed and simulated streamflow, are depicted in Figure 5.The range of the residuals is approximately [−13.7, 18.3] for ensemble runs and [−7.4,10.7] 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 Nash-Sutcliffe efficiency values are also indicated in the table.It is seen that the NS are very similar for all PSO modifications.The APartW variant achieved the smallest NS 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.minimum, where the error value less 10 −8 , in three functions.functions are  1 ,  2 (unimodal), and  7 (multimodal), and the results are comparable with our previous research [27].results were obtained by Liang and Suganthan [35], who applied distributed PSO algorithm with local search.They achieved the global in  1 ,  2  3 ,  6 , and  7 .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  1 ,  2 , and  4 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.

Discussion
Our results are similar to the results obtained by Nickabadi et al. [26].Their algorithms achieved the global optimum in functions  1 ,  2 , and  8 .
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  1 ,  2 , and  7 .Our results are better than the DE algorithm of Rönkkönen et al. [51], which achieved the global optimum in functions  1 and  7 or better than BLX MA method of Molina et al. [52], which converged only in functions  1 and  9 .

The Calibration of Bilan Model.
Our best obtained NS 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 rainfall-runoff model often used in sloping lowland catchments.The differences in their work are that, during the calibration of the model, they used 1 year of discharge observations and optimized 4 parameters.Their best achieved NS values for two catchments were 0.87 and 0.83.In our work, we used 30 years of observation and calibrated 6 parameters for 30 catchments.
Similar values of NS 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 234 catchments in the Czech Republic.The best achieved NS values during calibration of the Bilan model using SCDE were between 0.78 and 0.80, whereas using existing integrated local optimization technique with expert knowledge it was 0.73 [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.

Conclusions
The main aim of this paper was to test 4 selected PSO distributed versions on single-objective benchmark optimization problems and to apply them on calibration of hydrological model Bilan.For all 4 PSO versions, 3275 optimization problems were analysed, in which 275 minimizations for benchmark problems (i.e., 11 benchmark functions × 25 program runs) and 3000 inverse hydrological problems (i.e., 4 objective functions × 30 catchments × 25 program runs) were solved.
The new proposed variant APartW was compared with other existing PSO modifications-ConstrFactor, Lin-TimeVarW, and AdaptW on 11 benchmark functions prepared for the special session on real-parameter 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 ( 1 ,  2 , and  7 ).The APartW variant gave significantly better results in functions  4 ,  7 , and  11 .Beyond that, in functions  3 ,  5 , and  9 , there was no significant difference between the APartW and AdaptW methods.
All four PSO modifications were implemented into the Bilan rainfall-runoff 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 single-objective 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 rainfall-runoff simulations.

5. 1 . 4 :
The CEC 2005 Benchmark Optimization Problems.Both AdaptW and APartW modifications obtained the global Scatter plot of observed and runoff using APartW optimization.The optimized objective was NS.Catchment 01531000
[22,23]ticle'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 [− ⃗  max , ⃗  max ], where ⃗  max is equal to ⃗  max[22,23].

Table 2 :
lists Characteristics of the catchments.

Table 3 :
Statistical indices of 11 benchmark functions for the APartW and AdaptW modifications. 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 0.01 <  value ≤ 0.05; * * * for  value ≤ 0.001). 1

Table 4 :
Best values of all objective functions achieved in the total set of 30 catchments.The best calculated value for each optimized objective function (OOF) is highlighted in bold.

Table 5 :
The contrast test of the unadjusted medians with ranking.Rank is ranking based on contrast test; W.Rank is ranking based on Wilcoxon pair test.

Table 6 :
Parameters of Bilan model for catchment 01531000 obtained by the best model of each PSO modification.The optimized objective function was NS, and the last column indicates the value of that criterion.