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

Authors: Michala Jakubcová; Petr Máca; Pavel Pech
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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 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 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.


---

## Body

## 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 real-life 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 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 stand-alone 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. Section2 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.


## 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 velocityV
→
i
=
v
i
1
,
v
i
2
,
…
,
v
i
D
i
m is updated as (1)
V
→
i
t
+
1
=
V
→
i
t
+
c
1
·
U
→
1
t
⊗
P
→
i
t
-
X
→
i
t
+
c
2
·
U
→
2
t
⊗
G
→
t
-
X
→
i
t
,and the location X
→
i
=
x
i
1
,
x
i
2
,
…
,
x
i
D
i
m is simply defined as (2)
X
→
i
t
+
1
=
X
→
i
t
+
V
→
i
t
+
1
,for all i
=
1
,
…
,
S, where S is the total number of particles in swarm population, d
=
1
,
…
,
D
i
m, where D
i
m is the total number of dimensions (i.e., the number of parameters of hydrological model), c
1 and c
2 denote acceleration constants predefined by the user, and U
→
1 and U
→
2 are independent random vectors sampled from a uniform distribution in the range [
0,1
].The component with particle’s previous best positionP
→
i
=
p
i
1
,
p
i
2
,
…
,
p
i
D
i
m in (1) represents the cognition knowledge of swarm particles. When optimization problem is the minimization, f
(
P
→
i
t
)
≤
f
(
X
→
i
t
) for all so far known locations of given particle X
→
i
t. The component with the best position among all particles G
→
=
g
1
,
g
2
,
…
,
g
D
i
m in (1) controls the social influence of swarm particles. For all X
→
i in the population f
(
G
→
t
)
≤
f
(
X
→
i
t
). f is the analysed single-objective 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 [
-
V
→
m
a
x
,
V
→
m
a
x
], where V
→
m
a
x is equal to X
→
m
a
x [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 K is implemented into the PSO algorithm [24]. If the adapted particle’s velocity with K is as (3)
V
→
i
t
+
1
=
K
·
V
→
i
t
+
c
1
·
U
→
1
t
⊗
P
→
i
t
-
X
→
i
t
+
c
2
·
U
→
2
t
⊗
G
→
t
-
X
→
i
t
,then (4)
K
=
2
2
-
φ
-
φ
2
-
4
·
φ
,where φ
=
c
1
+
c
2 and φ
>
4.LinTimeVarW. In the second PSO modification, the linearly decreasing inertia weight w is used [25]. The adapted particle’s velocity equation is shown in (5)
V
→
i
t
+
1
=
w
·
V
→
i
t
+
c
1
·
U
→
1
t
⊗
P
→
i
t
-
X
→
i
t
+
c
2
·
U
→
2
t
⊗
G
→
t
-
X
→
i
t
.The value of inertia weight decreases at each iteration linearly from w
m
a
x
=
0.9 to w
m
i
n
=
0.4 and is simply defined as (6)
w
i
t
e
r
=
i
t
e
r
m
a
x
-
i
t
e
r
i
t
e
r
m
a
x
·
w
m
a
x
-
w
m
i
n
+
w
m
i
n
.AdaptW. One of the adaptive inertia weights from [26] was used in this paper. The modification of the inertia weight parameter is as (7)
w
i
t
e
r
=
w
m
a
x
-
w
m
i
n
·
P
s
→
t
+
w
m
i
n
,for which (8)
P
s
→
t
=
∑
i
=
1
n
S
→
i
t
n
,
S
→
i
t
=
1
if
f
P
→
i
t
<
f
P
→
i
t
-
1
0
if
f
P
→
i
t
=
f
P
→
i
t
-
1
,where S
→
i shows the success of ith particle, P
s
→ is the success percentage of the swarm, n is the size of the population, w
m
i
n
=
0, and w
m
a
x
=
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 f
(
X
→
i
t
)
<
f
(
P
→
i
t
-
1
), the local exploitation is supported. The inertia weight is calculated as (9)
w
→
i
t
=
w
m
a
x
+
w
m
i
n
2
-
w
m
i
n
·
U
→
t
+
w
m
i
n
,and thus, for random vector U
→ sampled from a uniform distribution in the range [
0,1
], w
m
i
n
=
0.1 and w
m
a
x
=
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, iff
(
X
→
i
t
)
≥
f
(
P
→
i
t
-
1
), the resulted inertia weight is computed as (10)
w
→
i
t
=
w
m
a
x
+
w
m
i
n
2
-
w
m
i
n
·
U
→
t
+
w
m
a
x
+
w
m
i
n
2
,where the variables are the same as in (9) and the inertia weight is thus between [
0.5,0.9
].


### 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 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 complexesN
C. The division is according to the functional value of each particle. All particles are sorted in increasing order, and then each jth complex receives X
→
j
,
X
→
j
+
N
C
,
X
→
j
+
2
N
C
,
… 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 of27 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 Algorithm1, where S
_
c
o
m
p is the number of particles in one complex. The other PSO modifications differ only on lines 9, 11, and 16 in Algorithm 1.Algorithm 1:PSO algorithm with APartW adaptation and SCE mechanism.
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 j
=
1  to  NC  do
(5)      A
j ← divided  E  into NC complexes
(6)      for  each particle i
=
1  to  S_comp in  A
j  do
(7)       if  (
f
(
X
→
i
)
<
f
(
P
→
i
)
)  then
(8)        P
→
i
=
X
→
i
(9)        w
→
i
=
w
m
a
x
+
w
m
i
n
2
-
w
m
i
n
⋅
U
→
+
w
m
i
n
(10)    else
(11)     w
→
i
=
w
m
a
x
+
w
m
i
n
2
-
w
m
i
n
⋅
U
→
+
w
m
a
x
+
w
m
i
n
2
(12)    end if
(13)   end for
(14)   for  each particle i
=
1  to  S_comp in  A
j  do
(15)    G
→
=
min
⁡
P
→
0
,
P
→
1
,
…
,
P
→
S
_
c
o
m
p
(16)    V
→
i
=
w
→
i
⋅
V
→
i
+
c
1
⋅
U
→
1
⊗
P
→
i
-
X
→
i
+
c
2
⋅
U
→
2
⊗
G
→
-
X
→
i
(17)    X
→
i
=
X
→
i
+
V
→
i
(18)   end for
(19)  end for
(20) until  termination criteria is met


## 2.1. The Original PSO

The original PSO algorithm estimates a new particle’s location using information of particle’s velocity. The velocityV
→
i
=
v
i
1
,
v
i
2
,
…
,
v
i
D
i
m is updated as (1)
V
→
i
t
+
1
=
V
→
i
t
+
c
1
·
U
→
1
t
⊗
P
→
i
t
-
X
→
i
t
+
c
2
·
U
→
2
t
⊗
G
→
t
-
X
→
i
t
,and the location X
→
i
=
x
i
1
,
x
i
2
,
…
,
x
i
D
i
m is simply defined as (2)
X
→
i
t
+
1
=
X
→
i
t
+
V
→
i
t
+
1
,for all i
=
1
,
…
,
S, where S is the total number of particles in swarm population, d
=
1
,
…
,
D
i
m, where D
i
m is the total number of dimensions (i.e., the number of parameters of hydrological model), c
1 and c
2 denote acceleration constants predefined by the user, and U
→
1 and U
→
2 are independent random vectors sampled from a uniform distribution in the range [
0,1
].The component with particle’s previous best positionP
→
i
=
p
i
1
,
p
i
2
,
…
,
p
i
D
i
m in (1) represents the cognition knowledge of swarm particles. When optimization problem is the minimization, f
(
P
→
i
t
)
≤
f
(
X
→
i
t
) for all so far known locations of given particle X
→
i
t. The component with the best position among all particles G
→
=
g
1
,
g
2
,
…
,
g
D
i
m in (1) controls the social influence of swarm particles. For all X
→
i in the population f
(
G
→
t
)
≤
f
(
X
→
i
t
). f is the analysed single-objective 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 [
-
V
→
m
a
x
,
V
→
m
a
x
], where V
→
m
a
x is equal to X
→
m
a
x [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 K is implemented into the PSO algorithm [24]. If the adapted particle’s velocity with K is as (3)
V
→
i
t
+
1
=
K
·
V
→
i
t
+
c
1
·
U
→
1
t
⊗
P
→
i
t
-
X
→
i
t
+
c
2
·
U
→
2
t
⊗
G
→
t
-
X
→
i
t
,then (4)
K
=
2
2
-
φ
-
φ
2
-
4
·
φ
,where φ
=
c
1
+
c
2 and φ
>
4.LinTimeVarW. In the second PSO modification, the linearly decreasing inertia weight w is used [25]. The adapted particle’s velocity equation is shown in (5)
V
→
i
t
+
1
=
w
·
V
→
i
t
+
c
1
·
U
→
1
t
⊗
P
→
i
t
-
X
→
i
t
+
c
2
·
U
→
2
t
⊗
G
→
t
-
X
→
i
t
.The value of inertia weight decreases at each iteration linearly from w
m
a
x
=
0.9 to w
m
i
n
=
0.4 and is simply defined as (6)
w
i
t
e
r
=
i
t
e
r
m
a
x
-
i
t
e
r
i
t
e
r
m
a
x
·
w
m
a
x
-
w
m
i
n
+
w
m
i
n
.AdaptW. One of the adaptive inertia weights from [26] was used in this paper. The modification of the inertia weight parameter is as (7)
w
i
t
e
r
=
w
m
a
x
-
w
m
i
n
·
P
s
→
t
+
w
m
i
n
,for which (8)
P
s
→
t
=
∑
i
=
1
n
S
→
i
t
n
,
S
→
i
t
=
1
if
f
P
→
i
t
<
f
P
→
i
t
-
1
0
if
f
P
→
i
t
=
f
P
→
i
t
-
1
,where S
→
i shows the success of ith particle, P
s
→ is the success percentage of the swarm, n is the size of the population, w
m
i
n
=
0, and w
m
a
x
=
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 f
(
X
→
i
t
)
<
f
(
P
→
i
t
-
1
), the local exploitation is supported. The inertia weight is calculated as (9)
w
→
i
t
=
w
m
a
x
+
w
m
i
n
2
-
w
m
i
n
·
U
→
t
+
w
m
i
n
,and thus, for random vector U
→ sampled from a uniform distribution in the range [
0,1
], w
m
i
n
=
0.1 and w
m
a
x
=
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, iff
(
X
→
i
t
)
≥
f
(
P
→
i
t
-
1
), the resulted inertia weight is computed as (10)
w
→
i
t
=
w
m
a
x
+
w
m
i
n
2
-
w
m
i
n
·
U
→
t
+
w
m
a
x
+
w
m
i
n
2
,where the variables are the same as in (9) and the inertia weight is thus between [
0.5,0.9
].


## 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 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 complexesN
C. The division is according to the functional value of each particle. All particles are sorted in increasing order, and then each jth complex receives X
→
j
,
X
→
j
+
N
C
,
X
→
j
+
2
N
C
,
… 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 of27 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 Algorithm1, where S
_
c
o
m
p is the number of particles in one complex. The other PSO modifications differ only on lines 9, 11, and 16 in Algorithm 1.Algorithm 1:PSO algorithm with APartW adaptation and SCE mechanism.
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 j
=
1  to  NC  do
(5)      A
j ← divided  E  into NC complexes
(6)      for  each particle i
=
1  to  S_comp in  A
j  do
(7)       if  (
f
(
X
→
i
)
<
f
(
P
→
i
)
)  then
(8)        P
→
i
=
X
→
i
(9)        w
→
i
=
w
m
a
x
+
w
m
i
n
2
-
w
m
i
n
⋅
U
→
+
w
m
i
n
(10)    else
(11)     w
→
i
=
w
m
a
x
+
w
m
i
n
2
-
w
m
i
n
⋅
U
→
+
w
m
a
x
+
w
m
i
n
2
(12)    end if
(13)   end for
(14)   for  each particle i
=
1  to  S_comp in  A
j  do
(15)    G
→
=
min
⁡
P
→
0
,
P
→
1
,
…
,
P
→
S
_
c
o
m
p
(16)    V
→
i
=
w
→
i
⋅
V
→
i
+
c
1
⋅
U
→
1
⊗
P
→
i
-
X
→
i
+
c
2
⋅
U
→
2
⊗
G
→
-
X
→
i
(17)    X
→
i
=
X
→
i
+
V
→
i
(18)   end for
(19)  end for
(20) until  termination criteria is met


## 3. 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 was15.The first set is represented by the11 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 of120 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].


### 3.1. The CEC 2005 Benchmark Optimization Problems

In order to compare the optimization algorithms, the single-objective 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 (f
1), Schwefel 1.2 (f
2), elliptic rotated (f
3), Schwefel 1.2 with noise (f
4), Schwefel 2.6 (f
5), Rosenbrock (f
6), Griewank rotated (f
7), Ackley rotated (f
8), Rastrigin (f
9), Rastrigin rotated (f
10), and Weierstrass rotated (f
11). The optimization problems f
1
–
f
5 are unimodal functions; the remaining f
6
–
f
11 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 5, and the maximum number of function evaluation is set to 10
000
·
D
i
m. 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 c
1
=
c
2
=
2, and the acceleration constants for adaptation using constriction factor are c
1
=
c
2
=
1.49445 with constriction coefficient K
=
0.729.


### 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 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 to3. 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.


#### 3.2.1. 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–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 first-order 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].Figure 1
Scheme of the Bilan rainfall-runoff model [18].The total amount of parameters of the daily version of the Bilan model is six, and they are listed in Table1. 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.Table 1
Parameters of the daily version of Bilan model [18].
Name
Description
PC1
Spa
Capacity of the soil moisture storage
Spa∈
〈
0,200
〉
Alf
Parameter of rainfall-runoff equation (direct runoff)
Alf∈
〈
0,1
〉
Dgm
Temperature/snow melting factor
Dgm∈
〈
0,200
〉
Soc
Parameter controlling distribution of percolation into interflow and groundwater recharge under summer conditions
Soc∈
〈
0,1
〉
Mec
Parameter controlling distribution of percolation into interflow and groundwater recharge under snow melt conditions
Mec∈
〈
0,1
〉
Grd
Parameter controlling outflow from groundwater storage (baseflow)
Grd∈
〈
0,1
〉
1PC = 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 Nash-Sutcliffe efficiency (NS) [32, 33, 46].The analysed objective functions are defined as(11)
MSE
=
1
N
∑
i
=
1
N
R
i
-
R
M
i
2
,
MAE
=
1
N
∑
i
=
1
N
R
i
-
R
M
i
,
MAPE
=
1
N
∑
i
=
1
N
R
i
-
R
M
i
R
i
,
FRV
=
∑
i
=
1
N
R
i
-
R
M
i
2
∑
i
=
1
N
R
i
-
R
¯
2
,where R is observed streamflow [mm], R
M is modelled streamflow [mm], R
¯ is mean of the observed streamflow [mm], and N is the total number of observations.The fraction residual variance (F
R
V) 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]. N
S is defined as (12)
N
S
=
1
-
F
R
V
.


#### 3.2.3. Dataset

For evaluation of the optimization ability of the proposed PSO modifications, we calibrated Bilan using datasets from30 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 from1948 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 lists the major characteristics of each catchment including latitude, longitude, drainage area, dominant soil, and vegetation types.Table 2
Characteristics of the catchments.
USGS ID
Lat.
Long.
Area[km2
]
Soil type
Veg. type1
01127000
41.5980
−71.9850
1147
Sandy loam
MF
01197500
42.2320
−73.3550
454
Sandy loam
MF
01321000
43.3528
−74.2708
790
Sandy loam
MF/CS
01371500
41.6860
−74.1660
1144
Silt loam
MF
01372500
41.6531
−73.8731
291
Silt loam
MF
01426500
42.0031
−75.3839
957
Silt loam
MF
01445500
40.8306
−74.9786
171
Sandy loam
MF
01503000
42.0353
−75.8033
3591
Silt loam
MF
01512500
42.2181
−75.8486
2386
Silt loam
MF
01518000
41.9083
−77.1297
454
Silt loam
MF
01531000
42.0022
−76.6350
4032
Silt loam
MF
01534000
41.5583
−75.8950
616
Loam
MF
01541000
40.8969
−78.6772
507
Loam
MF
01541500
40.9717
−78.4061
597
Loam
MF
01543500
41.3172
−78.1033
1102
Loam
MF
01548500
41.5217
−77.4478
972
Silt loam
MF
01556000
40.4631
−78.2000
468
Sandy loam
MF
01558000
40.6125
−78.1408
354
Sandy loam
MF
01559000
40.4850
−78.0190
1313
Sandy loam
MF
01560000
40.0717
−78.4928
277
Silt loam
MF
01562000
40.2158
−78.2656
1216
Silt loam
MF
01567000
40.4783
−77.1294
5397
Sandy loam
MF
01574000
40.0822
−76.7203
821
Silt loam
MF
01610000
39.5369
−78.4578
5002
Loam
MF
01628500
38.3220
−78.7550
1744
Clay
MF/CS
01631000
38.9139
−78.2111
2642
Loam
CS
01634000
38.9767
−78.3364
1236
Loam
CS
01643000
39.3880
−77.3800
1315
Silt loam
MF
01664000
38.5306
−77.8139
998
Clay loam
CS
01668000
38.3222
−77.5181
2568
Clay loam
CS
1MF = mixed forest; CS = closed shrublands.


## 3.1. The CEC 2005 Benchmark Optimization Problems

In order to compare the optimization algorithms, the single-objective 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 (f
1), Schwefel 1.2 (f
2), elliptic rotated (f
3), Schwefel 1.2 with noise (f
4), Schwefel 2.6 (f
5), Rosenbrock (f
6), Griewank rotated (f
7), Ackley rotated (f
8), Rastrigin (f
9), Rastrigin rotated (f
10), and Weierstrass rotated (f
11). The optimization problems f
1
–
f
5 are unimodal functions; the remaining f
6
–
f
11 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 5, and the maximum number of function evaluation is set to 10
000
·
D
i
m. 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 c
1
=
c
2
=
2, and the acceleration constants for adaptation using constriction factor are c
1
=
c
2
=
1.49445 with constriction coefficient K
=
0.729.


## 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 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 to3. 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.


### 3.2.1. 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–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 first-order 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].Figure 1
Scheme of the Bilan rainfall-runoff model [18].The total amount of parameters of the daily version of the Bilan model is six, and they are listed in Table1. 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.Table 1
Parameters of the daily version of Bilan model [18].
Name
Description
PC1
Spa
Capacity of the soil moisture storage
Spa∈
〈
0,200
〉
Alf
Parameter of rainfall-runoff equation (direct runoff)
Alf∈
〈
0,1
〉
Dgm
Temperature/snow melting factor
Dgm∈
〈
0,200
〉
Soc
Parameter controlling distribution of percolation into interflow and groundwater recharge under summer conditions
Soc∈
〈
0,1
〉
Mec
Parameter controlling distribution of percolation into interflow and groundwater recharge under snow melt conditions
Mec∈
〈
0,1
〉
Grd
Parameter controlling outflow from groundwater storage (baseflow)
Grd∈
〈
0,1
〉
1PC = 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 Nash-Sutcliffe efficiency (NS) [32, 33, 46].The analysed objective functions are defined as(11)
MSE
=
1
N
∑
i
=
1
N
R
i
-
R
M
i
2
,
MAE
=
1
N
∑
i
=
1
N
R
i
-
R
M
i
,
MAPE
=
1
N
∑
i
=
1
N
R
i
-
R
M
i
R
i
,
FRV
=
∑
i
=
1
N
R
i
-
R
M
i
2
∑
i
=
1
N
R
i
-
R
¯
2
,where R is observed streamflow [mm], R
M is modelled streamflow [mm], R
¯ is mean of the observed streamflow [mm], and N is the total number of observations.The fraction residual variance (F
R
V) 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]. N
S is defined as (12)
N
S
=
1
-
F
R
V
.


### 3.2.3. Dataset

For evaluation of the optimization ability of the proposed PSO modifications, we calibrated Bilan using datasets from30 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 from1948 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 lists the major characteristics of each catchment including latitude, longitude, drainage area, dominant soil, and vegetation types.Table 2
Characteristics of the catchments.
USGS ID
Lat.
Long.
Area[km2
]
Soil type
Veg. type1
01127000
41.5980
−71.9850
1147
Sandy loam
MF
01197500
42.2320
−73.3550
454
Sandy loam
MF
01321000
43.3528
−74.2708
790
Sandy loam
MF/CS
01371500
41.6860
−74.1660
1144
Silt loam
MF
01372500
41.6531
−73.8731
291
Silt loam
MF
01426500
42.0031
−75.3839
957
Silt loam
MF
01445500
40.8306
−74.9786
171
Sandy loam
MF
01503000
42.0353
−75.8033
3591
Silt loam
MF
01512500
42.2181
−75.8486
2386
Silt loam
MF
01518000
41.9083
−77.1297
454
Silt loam
MF
01531000
42.0022
−76.6350
4032
Silt loam
MF
01534000
41.5583
−75.8950
616
Loam
MF
01541000
40.8969
−78.6772
507
Loam
MF
01541500
40.9717
−78.4061
597
Loam
MF
01543500
41.3172
−78.1033
1102
Loam
MF
01548500
41.5217
−77.4478
972
Silt loam
MF
01556000
40.4631
−78.2000
468
Sandy loam
MF
01558000
40.6125
−78.1408
354
Sandy loam
MF
01559000
40.4850
−78.0190
1313
Sandy loam
MF
01560000
40.0717
−78.4928
277
Silt loam
MF
01562000
40.2158
−78.2656
1216
Silt loam
MF
01567000
40.4783
−77.1294
5397
Sandy loam
MF
01574000
40.0822
−76.7203
821
Silt loam
MF
01610000
39.5369
−78.4578
5002
Loam
MF
01628500
38.3220
−78.7550
1744
Clay
MF/CS
01631000
38.9139
−78.2111
2642
Loam
CS
01634000
38.9767
−78.3364
1236
Loam
CS
01643000
39.3880
−77.3800
1315
Silt loam
MF
01664000
38.5306
−77.8139
998
Clay loam
CS
01668000
38.3222
−77.5181
2568
Clay loam
CS
1MF = mixed forest; CS = closed shrublands.


## 3.2.1. 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–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 first-order 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].Figure 1
Scheme of the Bilan rainfall-runoff model [18].The total amount of parameters of the daily version of the Bilan model is six, and they are listed in Table1. 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.Table 1
Parameters of the daily version of Bilan model [18].
Name
Description
PC1
Spa
Capacity of the soil moisture storage
Spa∈
〈
0,200
〉
Alf
Parameter of rainfall-runoff equation (direct runoff)
Alf∈
〈
0,1
〉
Dgm
Temperature/snow melting factor
Dgm∈
〈
0,200
〉
Soc
Parameter controlling distribution of percolation into interflow and groundwater recharge under summer conditions
Soc∈
〈
0,1
〉
Mec
Parameter controlling distribution of percolation into interflow and groundwater recharge under snow melt conditions
Mec∈
〈
0,1
〉
Grd
Parameter controlling outflow from groundwater storage (baseflow)
Grd∈
〈
0,1
〉
1PC = 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 Nash-Sutcliffe efficiency (NS) [32, 33, 46].The analysed objective functions are defined as(11)
MSE
=
1
N
∑
i
=
1
N
R
i
-
R
M
i
2
,
MAE
=
1
N
∑
i
=
1
N
R
i
-
R
M
i
,
MAPE
=
1
N
∑
i
=
1
N
R
i
-
R
M
i
R
i
,
FRV
=
∑
i
=
1
N
R
i
-
R
M
i
2
∑
i
=
1
N
R
i
-
R
¯
2
,where R is observed streamflow [mm], R
M is modelled streamflow [mm], R
¯ is mean of the observed streamflow [mm], and N is the total number of observations.The fraction residual variance (F
R
V) 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]. N
S is defined as (12)
N
S
=
1
-
F
R
V
.


## 3.2.3. Dataset

For evaluation of the optimization ability of the proposed PSO modifications, we calibrated Bilan using datasets from30 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 from1948 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 lists the major characteristics of each catchment including latitude, longitude, drainage area, dominant soil, and vegetation types.Table 2
Characteristics of the catchments.
USGS ID
Lat.
Long.
Area[km2
]
Soil type
Veg. type1
01127000
41.5980
−71.9850
1147
Sandy loam
MF
01197500
42.2320
−73.3550
454
Sandy loam
MF
01321000
43.3528
−74.2708
790
Sandy loam
MF/CS
01371500
41.6860
−74.1660
1144
Silt loam
MF
01372500
41.6531
−73.8731
291
Silt loam
MF
01426500
42.0031
−75.3839
957
Silt loam
MF
01445500
40.8306
−74.9786
171
Sandy loam
MF
01503000
42.0353
−75.8033
3591
Silt loam
MF
01512500
42.2181
−75.8486
2386
Silt loam
MF
01518000
41.9083
−77.1297
454
Silt loam
MF
01531000
42.0022
−76.6350
4032
Silt loam
MF
01534000
41.5583
−75.8950
616
Loam
MF
01541000
40.8969
−78.6772
507
Loam
MF
01541500
40.9717
−78.4061
597
Loam
MF
01543500
41.3172
−78.1033
1102
Loam
MF
01548500
41.5217
−77.4478
972
Silt loam
MF
01556000
40.4631
−78.2000
468
Sandy loam
MF
01558000
40.6125
−78.1408
354
Sandy loam
MF
01559000
40.4850
−78.0190
1313
Sandy loam
MF
01560000
40.0717
−78.4928
277
Silt loam
MF
01562000
40.2158
−78.2656
1216
Silt loam
MF
01567000
40.4783
−77.1294
5397
Sandy loam
MF
01574000
40.0822
−76.7203
821
Silt loam
MF
01610000
39.5369
−78.4578
5002
Loam
MF
01628500
38.3220
−78.7550
1744
Clay
MF/CS
01631000
38.9139
−78.2111
2642
Loam
CS
01634000
38.9767
−78.3364
1236
Loam
CS
01643000
39.3880
−77.3800
1315
Silt loam
MF
01664000
38.5306
−77.8139
998
Clay loam
CS
01668000
38.3222
−77.5181
2568
Clay loam
CS
1MF = mixed forest; CS = closed shrublands.


## 4. Results

### 4.1. The CEC 2005 Benchmark Optimization Problems

Figure2 presents the convergence graphs for each benchmark function while utilizing all four distributed modifications of the PSO algorithm. The x-axis indicates the number of function evaluations and the y-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.Figure 2
Convergence graphs for the best fitness of the ConstrFactor, LinTimeVarW, AdaptW, and APartW modifications for the evaluated benchmark functionsf
1
–
f
11.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.Table3 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 P 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.Table 3
Statistical indices of11 benchmark functions for the APartW and AdaptW modifications.
Func.
PSO
Min
25%
Median
75%
Max
Mean
Std
P value1
f
1
APartW
5.68
E
-
14
3.98
E
-
13
6.93
E
-
12
8.41
E
-
10
1.94
E
-
05
8.52
E
-
07
3.87
E
-
06
1.000
AdaptW
0.00
E
+
00
5.68
E
-
14
1.14
E
-
13
1.71
E
-
13
2.27
E
-
13
1.02
E
-
13
6.14
E
-
14
f
2
APartW
5.68
E
-
14
5.12
E
-
13
9.66
E
-
13
5.46
E
-
12
2.68
E
-
11
3.73
E
-
12
5.77
E
-
12
0.970
AdaptW
5.68
E
-
14
2.27
E
-
13
3.98
E
-
13
6.82
E
-
13
1.53
E
-
12
5.16
E
-
13
4.33
E
-
13
f
3
APartW
3.72
E
+
05
1.99
E
+
06
4.05
E
+
06
7.30
E
+
06
2.92
E
+
07
6.33
E
+
06
6.71
E
+
06
0.439
AdaptW
4.17
E
+
05
7.08
E
+
05
1.02
E
+
06
1.57
E
+
06
2.30
E
+
06
1.15
E
+
06
5.78
E
+
05
f
4
APartW
1.48
E
-
01
3.29
E
+
00
9.75
E
+
00
5.00
E
+
01
1.17
E
+
03
1.02
E
+
02
2.48
E
+
02
0.000
∗
∗
∗
AdaptW
1.08
E
+
03
4.28
E
+
03
6.66
E
+
03
1.07
E
+
04
1.74
E
+
04
7.62
E
+
03
4.22
E
+
03
f
5
APartW
2.58
E
+
03
3.32
E
+
03
4.15
E
+
03
4.78
E
+
03
6.10
E
+
03
4.11
E
+
03
9.63
E
+
02
0.424
AdaptW
2.20
E
+
03
3.51
E
+
03
4.07
E
+
03
4.92
E
+
03
6.68
E
+
03
4.26
E
+
03
1.28
E
+
03
f
6
APartW
3.25
E
+
02
1.01
E
+
05
1.35
E
+
08
1.21
E
+
09
1.03
E
+
10
1.19
E
+
09
2.51
E
+
09
1.000
AdaptW
1.71
E
-
03
3.83
E
+
00
4.06
E
+
00
7.47
E
+
00
1.80
E
+
01
5.76
E
+
00
4.87
E
+
00
f
7
APartW
2.84
E
-
14
3.69
E
-
11
1.23
E
-
02
3.45
E
-
02
7.72
E
+
03
3.09
E
+
02
1.54
E
+
03
0.040
∗
AdaptW
2.84
E
-
14
9.86
E
-
03
2.22
E
-
02
3.69
E
-
02
6.15
E
-
02
2.56
E
-
02
1.69
E
-
02
f
8
APartW
2.08
E
+
01
2.16
E
+
01
2.17
E
+
01
2.18
E
+
01
2.19
E
+
01
2.17
E
+
01
2.03
E
-
01
0.968
AdaptW
2.06
E
+
01
2.14
E
+
01
2.16
E
+
01
2.18
E
+
01
2.20
E
+
01
2.15
E
+
01
3.89
E
-
01
f
9
APartW
6.37
E
+
01
8.66
E
+
01
1.08
E
+
02
1.21
E
+
02
1.56
E
+
02
1.05
E
+
02
2.58
E
+
01
0.770
AdaptW
5.67
E
+
01
7.36
E
+
01
9.35
E
+
01
1.17
E
+
02
1.50
E
+
02
9.48
E
+
01
2.57
E
+
01
f
10
APartW
1.60
E
+
02
7.92
E
+
02
1.23
E
+
03
1.41
E
+
03
1.96
E
+
03
1.14
E
+
03
5.13
E
+
02
0.973
AdaptW
1.34
E
+
02
2.05
E
+
02
2.46
E
+
02
3.03
E
+
02
6.71
E
+
02
2.59
E
+
02
1.05
E
+
02
f
11
APartW
2.73
E
+
01
3.16
E
+
01
3.33
E
+
01
3.58
E
+
01
3.94
E
+
01
3.35
E
+
01
3.30
E
+
00
0.049
∗
AdaptW
2.93
E
+
01
3.44
E
+
01
3.72
E
+
01
3.99
E
+
01
6.30
E
+
01
4.13
E
+
01
1.08
E
+
01
1
P 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
<
P value ≤ 0.05; ∗
∗
∗ for P value ≤ 0.001).The results of the analysis in Table3 show that the APartW modification gives significantly better results for three benchmark functions (f
4, f
7, and f
11). In functions f
3, f
5, and f
9, there is no significant difference between APartW and AdaptW. Beyond that, in functions f
1 and f
2, both APartW and AdaptW found the global minimum. For multimodal functions f
6, f
8, and f
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.


### 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 ID01531000) 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 to25. 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 Table4 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.Table 4
Best values of all objective functions achieved in the total set of30 catchments. The best calculated value for each optimized objective function (OOF) is highlighted in bold.
OOF
PSO modif.
MSE
MAE
MAPE
NS
MSE
ConstrFactor
6.674
E
-
01
3.473
E
-
01
4.186
E
-
01
7.172
E
-
01
LinTimeVarW
5.816
E
-
01
3.598
E
-
01
4.041
E
-
01
7.466
E
-
01
AdaptW
5.816
E
-
01
3.594
E
-
01
4.144
E
-
01
7.466
E
-
01
APartW
6.253
E
-
01
3.547
E
-
01
3.376
E
-
01
7.297
E
-
01
MAE
ConstrFactor
7.432
E
-
01
3.564
E
-
01
3.380
E
-
01
7.097
E
-
01
LinTimeVarW
7.776
E
-
01
3.445
E
-
01
3.380
E
-
01
7.339
E
-
01
AdaptW
6.822
E
-
01
3.421
E
-
01
3.420
E
-
01
7.329
E
-
01
APartW
6.854
E
-
01
3.442
E
-
01
3.484
E
-
01
7.259
E
-
01
MAPE
ConstrFactor
7.800
E
-
01
3.550
E
-
01
3.290
E
-
01
6.808
E
-
01
LinTimeVarW
8.020
E
-
01
3.477
E
-
01
3.181
E
-
01
6.962
E
-
01
AdaptW
8.020
E
-
01
3.480
E
-
01
3.238
E
-
01
6.894
E
-
01
APartW
6.849
E
-
01
3.485
E
-
01
3.300
E
-
01
7.202
E
-
01
NS
ConstrFactor
7.268
E
-
01
3.544
E
-
01
3.380
E
-
01
7.257
E
-
01
LinTimeVarW
5.816
E
-
01
3.605
E
-
01
4.015
E
-
01
7.466
E
-
01
AdaptW
5.816
E
-
01
3.594
E
-
01
4.144
E
-
01
7.466
E
-
01
APartW
6.219
E
-
01
3.504
E
-
01
3.581
E
-
01
7.333
E
-
01In 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:M
S
E
=
1.044
E
+
00, M
A
E
=
3.598
E
-
01, M
A
P
E
=
3.287
E
-
01, and N
S
=
7.462
E
-
01. The results are for optimization based on M
S
E, M
A
E, M
A
P
E, and N
S, respectively. In comparison with Table 4, the optimization with PSO method is always better than with the integrated BS method.Table5 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.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.
ConstrFactor
LinTimeVarW
AdaptW
APartW
Rank
W.Rank
MSE
ConstrFactor
—
372.66
373.23
373.22
4
4
LinTimeVarW
−372.66
—
0.57
0.55
3
3
AdaptW
−373.23
−0.57
—
−0.01
1
2
APartW
−373.22
−0.55
0.01
—
2
1
MAE
ConstrFactor
—
374.78
374.82
374.82
4
4
LinTimeVarW
−374.78
—
0.04
0.04
3
3
AdaptW
−374.82
−0.04
—
0.00
2
2
APartW
−374.82
−0.04
−0.00
—
1
1
MAPE
ConstrFactor
—
374.96
375.02
375.02
4
4
LinTimeVarW
−374.96
—
0.06
0.05
3
3
AdaptW
−375.02
−0.06
—
−0.00
1
1-2
APartW
−375.02
−0.05
0.00
—
2
1-2
NS
ConstrFactor
—
375.01
374.91
374.91
4
4
LinTimeVarW
−375.01
—
−0.10
−0.10
1
1
AdaptW
−374.91
0.10
—
−0.00
2
2
APartW
−374.91
0.10
0.00
—
3
3In 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 M
S
E and M
A
P
E 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 Figure3 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 N
S.Figure 3
Observed streamflow and corresponding simulations from Bilan model using APartW optimization. The optimized objective function wasN
S. Catchment 01531000, year 1976.From Figure3 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.Figure 4
Scatter plot of observed and modelled runoff using APartW optimization. The optimized objective function wasN
S. Catchment 01531000, year 1976.Figure 5
Calculated residuals using APartW optimization. The optimized objective function wasN
S. Catchment 01531000, year 1976.The values of parameters obtained by the best model of each PSO modification for the chosen catchment are in Table6. The obtained Nash-Sutcliffe efficiency values are also indicated in the table. It is seen that the N
S  are very similar for all PSO modifications. The APartW variant achieved the smallest N
S value, but it was not significantly smaller than the others.Table 6
Parameters of Bilan model for catchment01531000 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.
PSO modif.
Spa
Alf
Dgm
Soc
Mec
Grd
NS
ConstrFactor
198.7
0.091
1.55
7.05
E
-
3
10.68
E
-
3
0.352
0.622
LinTimeVarW
199.9
0.559
1.59
2.33
E
-
3
6.00
E
-
3
0.276
0.644
AdaptW
200.0
0.552
1.59
2.38
E
-
3
6.09
E
-
3
0.275
0.644
APartW
178.3
0.186
1.33
5.85
E
-
3
17.33
E
-
3
0.301
0.621The 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.


## 4.1. The CEC 2005 Benchmark Optimization Problems

Figure2 presents the convergence graphs for each benchmark function while utilizing all four distributed modifications of the PSO algorithm. The x-axis indicates the number of function evaluations and the y-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.Figure 2
Convergence graphs for the best fitness of the ConstrFactor, LinTimeVarW, AdaptW, and APartW modifications for the evaluated benchmark functionsf
1
–
f
11.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.Table3 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 P 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.Table 3
Statistical indices of11 benchmark functions for the APartW and AdaptW modifications.
Func.
PSO
Min
25%
Median
75%
Max
Mean
Std
P value1
f
1
APartW
5.68
E
-
14
3.98
E
-
13
6.93
E
-
12
8.41
E
-
10
1.94
E
-
05
8.52
E
-
07
3.87
E
-
06
1.000
AdaptW
0.00
E
+
00
5.68
E
-
14
1.14
E
-
13
1.71
E
-
13
2.27
E
-
13
1.02
E
-
13
6.14
E
-
14
f
2
APartW
5.68
E
-
14
5.12
E
-
13
9.66
E
-
13
5.46
E
-
12
2.68
E
-
11
3.73
E
-
12
5.77
E
-
12
0.970
AdaptW
5.68
E
-
14
2.27
E
-
13
3.98
E
-
13
6.82
E
-
13
1.53
E
-
12
5.16
E
-
13
4.33
E
-
13
f
3
APartW
3.72
E
+
05
1.99
E
+
06
4.05
E
+
06
7.30
E
+
06
2.92
E
+
07
6.33
E
+
06
6.71
E
+
06
0.439
AdaptW
4.17
E
+
05
7.08
E
+
05
1.02
E
+
06
1.57
E
+
06
2.30
E
+
06
1.15
E
+
06
5.78
E
+
05
f
4
APartW
1.48
E
-
01
3.29
E
+
00
9.75
E
+
00
5.00
E
+
01
1.17
E
+
03
1.02
E
+
02
2.48
E
+
02
0.000
∗
∗
∗
AdaptW
1.08
E
+
03
4.28
E
+
03
6.66
E
+
03
1.07
E
+
04
1.74
E
+
04
7.62
E
+
03
4.22
E
+
03
f
5
APartW
2.58
E
+
03
3.32
E
+
03
4.15
E
+
03
4.78
E
+
03
6.10
E
+
03
4.11
E
+
03
9.63
E
+
02
0.424
AdaptW
2.20
E
+
03
3.51
E
+
03
4.07
E
+
03
4.92
E
+
03
6.68
E
+
03
4.26
E
+
03
1.28
E
+
03
f
6
APartW
3.25
E
+
02
1.01
E
+
05
1.35
E
+
08
1.21
E
+
09
1.03
E
+
10
1.19
E
+
09
2.51
E
+
09
1.000
AdaptW
1.71
E
-
03
3.83
E
+
00
4.06
E
+
00
7.47
E
+
00
1.80
E
+
01
5.76
E
+
00
4.87
E
+
00
f
7
APartW
2.84
E
-
14
3.69
E
-
11
1.23
E
-
02
3.45
E
-
02
7.72
E
+
03
3.09
E
+
02
1.54
E
+
03
0.040
∗
AdaptW
2.84
E
-
14
9.86
E
-
03
2.22
E
-
02
3.69
E
-
02
6.15
E
-
02
2.56
E
-
02
1.69
E
-
02
f
8
APartW
2.08
E
+
01
2.16
E
+
01
2.17
E
+
01
2.18
E
+
01
2.19
E
+
01
2.17
E
+
01
2.03
E
-
01
0.968
AdaptW
2.06
E
+
01
2.14
E
+
01
2.16
E
+
01
2.18
E
+
01
2.20
E
+
01
2.15
E
+
01
3.89
E
-
01
f
9
APartW
6.37
E
+
01
8.66
E
+
01
1.08
E
+
02
1.21
E
+
02
1.56
E
+
02
1.05
E
+
02
2.58
E
+
01
0.770
AdaptW
5.67
E
+
01
7.36
E
+
01
9.35
E
+
01
1.17
E
+
02
1.50
E
+
02
9.48
E
+
01
2.57
E
+
01
f
10
APartW
1.60
E
+
02
7.92
E
+
02
1.23
E
+
03
1.41
E
+
03
1.96
E
+
03
1.14
E
+
03
5.13
E
+
02
0.973
AdaptW
1.34
E
+
02
2.05
E
+
02
2.46
E
+
02
3.03
E
+
02
6.71
E
+
02
2.59
E
+
02
1.05
E
+
02
f
11
APartW
2.73
E
+
01
3.16
E
+
01
3.33
E
+
01
3.58
E
+
01
3.94
E
+
01
3.35
E
+
01
3.30
E
+
00
0.049
∗
AdaptW
2.93
E
+
01
3.44
E
+
01
3.72
E
+
01
3.99
E
+
01
6.30
E
+
01
4.13
E
+
01
1.08
E
+
01
1
P 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
<
P value ≤ 0.05; ∗
∗
∗ for P value ≤ 0.001).The results of the analysis in Table3 show that the APartW modification gives significantly better results for three benchmark functions (f
4, f
7, and f
11). In functions f
3, f
5, and f
9, there is no significant difference between APartW and AdaptW. Beyond that, in functions f
1 and f
2, both APartW and AdaptW found the global minimum. For multimodal functions f
6, f
8, and f
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.


## 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 ID01531000) 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 to25. 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 Table4 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.Table 4
Best values of all objective functions achieved in the total set of30 catchments. The best calculated value for each optimized objective function (OOF) is highlighted in bold.
OOF
PSO modif.
MSE
MAE
MAPE
NS
MSE
ConstrFactor
6.674
E
-
01
3.473
E
-
01
4.186
E
-
01
7.172
E
-
01
LinTimeVarW
5.816
E
-
01
3.598
E
-
01
4.041
E
-
01
7.466
E
-
01
AdaptW
5.816
E
-
01
3.594
E
-
01
4.144
E
-
01
7.466
E
-
01
APartW
6.253
E
-
01
3.547
E
-
01
3.376
E
-
01
7.297
E
-
01
MAE
ConstrFactor
7.432
E
-
01
3.564
E
-
01
3.380
E
-
01
7.097
E
-
01
LinTimeVarW
7.776
E
-
01
3.445
E
-
01
3.380
E
-
01
7.339
E
-
01
AdaptW
6.822
E
-
01
3.421
E
-
01
3.420
E
-
01
7.329
E
-
01
APartW
6.854
E
-
01
3.442
E
-
01
3.484
E
-
01
7.259
E
-
01
MAPE
ConstrFactor
7.800
E
-
01
3.550
E
-
01
3.290
E
-
01
6.808
E
-
01
LinTimeVarW
8.020
E
-
01
3.477
E
-
01
3.181
E
-
01
6.962
E
-
01
AdaptW
8.020
E
-
01
3.480
E
-
01
3.238
E
-
01
6.894
E
-
01
APartW
6.849
E
-
01
3.485
E
-
01
3.300
E
-
01
7.202
E
-
01
NS
ConstrFactor
7.268
E
-
01
3.544
E
-
01
3.380
E
-
01
7.257
E
-
01
LinTimeVarW
5.816
E
-
01
3.605
E
-
01
4.015
E
-
01
7.466
E
-
01
AdaptW
5.816
E
-
01
3.594
E
-
01
4.144
E
-
01
7.466
E
-
01
APartW
6.219
E
-
01
3.504
E
-
01
3.581
E
-
01
7.333
E
-
01In 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:M
S
E
=
1.044
E
+
00, M
A
E
=
3.598
E
-
01, M
A
P
E
=
3.287
E
-
01, and N
S
=
7.462
E
-
01. The results are for optimization based on M
S
E, M
A
E, M
A
P
E, and N
S, respectively. In comparison with Table 4, the optimization with PSO method is always better than with the integrated BS method.Table5 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.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.
ConstrFactor
LinTimeVarW
AdaptW
APartW
Rank
W.Rank
MSE
ConstrFactor
—
372.66
373.23
373.22
4
4
LinTimeVarW
−372.66
—
0.57
0.55
3
3
AdaptW
−373.23
−0.57
—
−0.01
1
2
APartW
−373.22
−0.55
0.01
—
2
1
MAE
ConstrFactor
—
374.78
374.82
374.82
4
4
LinTimeVarW
−374.78
—
0.04
0.04
3
3
AdaptW
−374.82
−0.04
—
0.00
2
2
APartW
−374.82
−0.04
−0.00
—
1
1
MAPE
ConstrFactor
—
374.96
375.02
375.02
4
4
LinTimeVarW
−374.96
—
0.06
0.05
3
3
AdaptW
−375.02
−0.06
—
−0.00
1
1-2
APartW
−375.02
−0.05
0.00
—
2
1-2
NS
ConstrFactor
—
375.01
374.91
374.91
4
4
LinTimeVarW
−375.01
—
−0.10
−0.10
1
1
AdaptW
−374.91
0.10
—
−0.00
2
2
APartW
−374.91
0.10
0.00
—
3
3In 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 M
S
E and M
A
P
E 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 Figure3 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 N
S.Figure 3
Observed streamflow and corresponding simulations from Bilan model using APartW optimization. The optimized objective function wasN
S. Catchment 01531000, year 1976.From Figure3 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.Figure 4
Scatter plot of observed and modelled runoff using APartW optimization. The optimized objective function wasN
S. Catchment 01531000, year 1976.Figure 5
Calculated residuals using APartW optimization. The optimized objective function wasN
S. Catchment 01531000, year 1976.The values of parameters obtained by the best model of each PSO modification for the chosen catchment are in Table6. The obtained Nash-Sutcliffe efficiency values are also indicated in the table. It is seen that the N
S  are very similar for all PSO modifications. The APartW variant achieved the smallest N
S value, but it was not significantly smaller than the others.Table 6
Parameters of Bilan model for catchment01531000 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.
PSO modif.
Spa
Alf
Dgm
Soc
Mec
Grd
NS
ConstrFactor
198.7
0.091
1.55
7.05
E
-
3
10.68
E
-
3
0.352
0.622
LinTimeVarW
199.9
0.559
1.59
2.33
E
-
3
6.00
E
-
3
0.276
0.644
AdaptW
200.0
0.552
1.59
2.38
E
-
3
6.09
E
-
3
0.275
0.644
APartW
178.3
0.186
1.33
5.85
E
-
3
17.33
E
-
3
0.301
0.621The 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 than1
0
-
8, in three functions. These functions are f
1, f
2 (unimodal), and f
7 (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 f
1, f
2, f
3, f
6, and f
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 f
1, f
2, and f
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.Our results are similar to the results obtained by Nickabadi et al. [26]. Their algorithms achieved the global optimum in functions f
1, f
2, and f
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 f
1, f
2, and f
7. Our results are better than the DE algorithm of Rönkkönen et al. [51], which achieved the global optimum in functions f
1 and f
7 or better than BLX_MA method of Molina et al. [52], which converged only in functions f
1 and f
9.


### 5.2. The Calibration of Bilan Model

Our best obtainedN
S 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 N
S 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 ofN
S 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 on234 catchments in the Czech Republic. The best achieved N
S 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.


## 5.1. The CEC 2005 Benchmark Optimization Problems

Both AdaptW and APartW modifications obtained the global minimum, where the error value was less than1
0
-
8, in three functions. These functions are f
1, f
2 (unimodal), and f
7 (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 f
1, f
2, f
3, f
6, and f
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 f
1, f
2, and f
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.Our results are similar to the results obtained by Nickabadi et al. [26]. Their algorithms achieved the global optimum in functions f
1, f
2, and f
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 f
1, f
2, and f
7. Our results are better than the DE algorithm of Rönkkönen et al. [51], which achieved the global optimum in functions f
1 and f
7 or better than BLX_MA method of Molina et al. [52], which converged only in functions f
1 and f
9.


## 5.2. The Calibration of Bilan Model

Our best obtainedN
S 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 N
S 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 ofN
S 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 on234 catchments in the Czech Republic. The best achieved N
S 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.


## 6. Conclusions

The main aim of this paper was to test4 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, LinTimeVarW, and AdaptW on11 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 (f
1, f
2, and f
7). The APartW variant gave significantly better results in functions f
4, f
7, and f
11. Beyond that, in functions f
3, f
5, and f
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.



---

*Source: 968067-2015-10-04.xml*