Hybrid Cuckoo Search and Harmony Search Algorithm and Its

Modifications for the Calibration of Groundwater Flow Models

D. K. Valetov

, G. D. Neuvazhaev

, V. S. Svitelman

and E. A. Saveleva

Nuclear Safety Institute of the Russian Academy of Sciences, Bolshaya Tulskaya Street, Moscow, Russia

Keywords: Groundwater Flow Modelling, Model Calibration, Optimization Methods, Safety Assessment, Uncertainties.

Abstract: Due to inherent uncertainties associated with the groundwater system characterization, the model calibration

is the significant step for the obtaining of the reliable predictions that could be used in long-term safety

assessment. The following paper is focused on the modelling of the groundwater flow in heterogeneous media

using the data of the geological engineering survey for the prospective site of the radioactive waste deep

geological disposal at Nizhnekansky granite-gneiss crystalline rock massif (Krasnoyarsk Territory). This case

study illustrates the efficiency of heuristic optimization methods and hybridization approach for the

calibration of groundwater models.

1 INTRODUCTION

Safety assessments of the prospective radioactive

waste disposal are based on models of the disposal

facility and of its natural surroundings. Numerical

models of radionuclide transport processes in

geological media require the set of site-specific

parameters like flow and transport properties,

boundary conditions and so on. And these parameters

are associated with significant uncertainties due to the

natural variability of the geological media, lack of the

ability to measure them at each point of interest, and

the simplifications during the construction of the

numerical model.

Model calibration aims to reduce uncertainty in

the parametrization of the numerical model by

comparing the model predictions with site-specific

field observations and measurements. In practice, if a

model can be calibrated successfully for a variety of

site-specific conditions, it means an increased level of

confidence in the model’s ability to represent the

system behaviour and estimate its effects that cannot

be measured.

The model calibration has always been an

essential step for groundwater flow and transport

https://orcid.org/0000-0001-5432-4084

https://orcid.org/0000-0001-6821-1638

https://orcid.org/0000-0002-5976-6049

https://orcid.org/0000-0002-6562-8750

modelling, and the approaches and techniques have

been developing, evolving and inter-influencing:

from the manual calibration through gradient-based

and direct search to the heuristic and various hybrid

methods. Gradient-based methods are based on the

use of local derivative information to guide the search

direction and may not achieve an optimal solution in

complex problems with a large number of decision

variables due to converging to local suboptimal

solutions.

The alternatives are to utilize deterministic

pattern-based procedures (direct search) or adopt

structured randomness elements from natural

optimization strategies (heuristic and meta-heuristic

algorithms) (Maier et al., 2014). The last ones have

been extensively used lately to solve calibration

problems in different scientific fields including

groundwater flow and transport modelling because

they appeared to be able to overcome such common

challenges as mixed parameter types, the

nonlinearity, the discontinuities and the local minima

(Maier et al., 2014), (Haddad et al., 2013).

The problem under consideration in this paper is

parametrization of the cross-sectional groundwater

flow model in the heterogeneous geological media

Valetov, D., Neuvazhaev, G., Svitelman, V. and Saveleva, E.

Hybrid Cuckoo Search and Harmony Search Algorithm and Its Modiﬁcations for the Calibration of Groundwater Flow Models.

DOI: 10.5220/0008345502210228

In Proceedings of the 11th International Joint Conference on Computational Intelligence (IJCCI 2019), pages 221-228

ISBN: 978-989-758-384-1

221

Figure 1: Schematic geological cross-section of the modelling area, model 1.

Figure 2: Schematic geological cross-section of the modelling area, model 2.

based on the data of geological survey in the area of

the prospective radioactive waste deep geological

disposal. The parameter optimization tool used for

this model is a hybrid algorithm (Fister, 2015) that

combines the advantages of two different heuristic

methods. First, Cuckoo search (CS) has great

exploratory qualities due to the use of Levy Flight.

And the second one, Harmony search (HS), provides

efficient exploitation of the search space by efficient

adjustment of positions of the search agents. In

addition to the existent hybrid, two improvements of

it are considered: the use of CS with a chaotic varying

step size (Wang et al., 2016) and the modification of

the distribution for the HS.

2 PROBLEM DESIGN

Russian deep geological repository project has been

started about 30 years ago. To date, the siting

procedure is completed, and the license for the

construction of the underground research facility

(URF) at Nizhnekansky granite-gneiss crystalline

rock massif (Krasnoyarsk Territory) is received. The

scope and purposes of the URF are site

characterization, technology development and

testing, safety assessment and demonstration in

support of the development of the deep geological

disposal. The current concept under consideration of

the deep geological disposal is following: two

disposal levels will be located at 450-525 m depth.

High-level waste containers will be located in deep

vertical boreholes. That is consistent with

internationally considered mined deep borehole

matrix concept (Baldwin, Chapman and Neall, 2008).

During the geological survey of the site the wells

were drilled to depths of 700 m, the packer tests with

interval pressure measurements and single pumping

tests were carried out. The groundwater flow model

represents the steady-state problem in the confined

media.

The differential equations that describe ground-

ECTA 2019 - 11th International Conference on Evolutionary Computation Theory and Applications

222

water flow are the following:

Qu 

(1)

u K h  

(2)

where

is the velocity vector determined by Darcy

law (2) and h – hydraulic head, Q – sources and

drains, K – hydraulic conductivity.

These equations are solved by the finite-volume

method with two-point flow approximation by means

of GeRa groundwater flow and transport modelling

computer code (Kapyrin et al., 2018).

The two models have been constructed on the

basis of the structural geological cross-sections

corresponding to the two alternative interpretations of

the well data (Jobmann, 2016) (figures 1 and 2).

Different geological structural elements (a1 – a3 for

model 1 and b1 – b7 for model 2) are indicated by

different colors.

The boundary conditions for each model were set

in the following way:

 The right border was drawn along streamlines

ending on the Shumikha River (the streamlines

from the hydrodynamic point of view are the

impermeable boundaries);

 The lower boundary corresponds to no-flow

boundary condition due to the fact that the rocks

that lie at the mark -350m are treated as

impermeable;

 On the upper boundary the rainfall recharge rate

was set to 0.00001 m/day, except for the

Shumikha River on which Dirichlet-type

boundary condition with the level of 331m was

set;

 At the left border, Dirichlet-type boundary

condition with the constant hydraulic head of

400m was set, which was obtained on the basis of

water level measurements.

The hydraulic conductivity coefficients were varied

during calibration in the separate range for each

structural element (3 parameters for the model 1 and

7 parameters for model 2). Parameter variation ranges

were set by expert assessment and given in the

Results section. The underlying rationale is the

following. For intact rock specimen ranges of the

hydraulic conductivities were set according to the

geological survey findings (Jobmann, 2016). And for

fractured rock (b3.Fissured dikes, b5.Fissured gneiss,

b6.Сrushing zone) were set wider ranges because

fractures could be both water-conductive and non-

conductive based on the available information.

Observations for the calibration were represented

by 35 known hydraulic heads from 3 pumping wells.

3 CALIBRATION AS

OPTIMIZATION PROBLEM

For credible safety assessment of the geological

disposal system a model that properly represents an

actual system is required. Hence the model

calibration procedure aims to obtain a set of

parameters that would produce the best possible fit of

simulated and available observed values quantified

by the target function (Hill, 2004). During this

procedure, the vector of model input parameters is

varied by a specific algorithm, and the target function

(most commonly it is sum of weighted least-squares)

is evaluated at each step.

Figure 3: Scheme of the interaction between black box

model and optimization algorithm.

Therefore, the calibration is a particular case of

the global optimization problem where groundwater

flow is represented as a black box with the vector of

the input parameter and the target function as the

response variable (figure 3). The aim is to find the

minimum of target function surface in parametric

search space.

Such characteristics of groundwater flow and

transport models as nonlinearity, multimodality, the

high dimensionality of search space makes the

application of the gradient-based methods inefficient

for them. These conditions made the use of the

heuristic algorithms that combine rules and

randomness to imitate natural phenomena very

popular. The reason for the variety of heuristic

algorithms is explained by «no-free-lunch theorem of

optimization» that states that a general-purpose,

universal optimization strategy is impossible.

The only way one heuristic method can

outperform another is if it is specialized to the

structure of the specific problem. The promising way

of such kind of specialization is hybrids that combine

multiple algorithms in order to employ their

advantages and compensate for the weaknesses. It is

Hybrid Cuckoo Search and Harmony Search Algorithm and Its Modiﬁcations for the Calibration of Groundwater Flow Models

223

Table 1: Adjusted properties of the optimization methods.

Algorithm

Hyper-parameter

Value

Cuckoo Search

Population size (CMS)

Abandon probability (PA)

0.1

Step size scaling factor (α)

0.25

Dispersion coefficient for Levy Flight

1.5

Harmony Search

Population size (HMS)

Harmony Memory Considering Rate (HMCR)

0.5

Harmony Memory Considering Rate (HMCR) (in hybrid

algorithm)

1.0

Pitch adjusting rate (PAR)

0.5

Step size scaling factor (α)

0.01

Hybrid parameters

Number of Harmony Search iterations to be done per each

Cuckoo Search iteration

Chaos map to be used in Chaos Cuckoo Search

logistic

Dispersion coefficient for distribution used in modified

Harmony Search

6.67

worth mentioning that the concept of hybrids can be

found in the engineering of nature-based solutions as

is the case, for example, in the human brain where

applying multiple memory systems guarantees the

optimality of decision-making (Hamid and Braun,

2019).

The hybrid algorithm (CS/HS) that combines

Cuckoo Search and Harmony Search was proposed

by Wang et al. (Wang et al., 2016). The basic

algorithm in the hybrid is CS (Gandomi, Yang and

Alavi, 2013). The CS algorithm is inspired by the

parasitic brood reproductive strategy of cuckoo birds.

Specifics of cuckoos nesting process are following:

cuckoos lay their eggs in the nests of other birds; the

host bird could discover them with some probability

(PA) and abandon the nest or throw the alien egg

away. In the algorithm each cuckoo is represented by

a point in search space, if an egg has not been

discovered by the host bird, then it remains and

another cuckoo appears at this point. Target function

in this metaphor is the quality of the nest (lack of

predators, food availability, etc.)

The CS algorithm uses Levy flight search pattern

to explore search space. It mimics the behavior that is

very common for different bird and insects species:

series of straight flight paths punctuated by sudden

turns (Cole, 1995), (Rhee et al., 2011). Such

characteristics of Levy distribution as infinite

mathematical expectation and variation, and power-

law tail and positive asymmetry have a positive effect

on global search in the search space. Search process

with step sampled from this distribution covers the

search space faster than if with normal distribution

used. Moreover, bell-bottomness of this distribution

advances sampling of medium steps and thus better

local search too. This optimization method is

effectively used for high dimension problems in a

wide range of problem fields, including geological

modelling problems, for example, underground water

resources estimate (Gupta, Das and Panchal, 2013).

To sum up, the hyper-parameters (the options that

can customize the algorithm’s behavior) are the

following:

 population size (Cuckoo Memory Size – CMS);

 probability for an egg to be discovered by host

(Abandon Probability – PA);

 step size scaling factor related to the scales of the

problem of interest (α);

 dispersion coefficient for Levy Flight distribution

(σ) (Cole, 1995), (Rhee et al., 2011);

 step size scaling factor related to the scales of the

problem of interest (α).

The second component of the hybrid is HS (Ayvaz,

2009). It reproduces the idealized natural musical

improvisation processes. In this metaphor parameter

values are pitches of instruments, the target function

is the aesthetic standard, and the global optimum is

perfect harmony.

The search process includes three components:

memory consideration, pitch adjustment, and random

selection. Properties of the algorithm are the

following:

 amount of harmonies (parameter vectors) in

memory (Harmony Memory Size – HMS);

 probability to choose an instrument pitch from

memory (Harmony Memory Considering Rate –

HMCR);

 the probability to adjust harmony chosen from

memory (Pitch Adjusting Rate – PAR);

 step size scaling factor related to the scales of the

problem of interest (α).

ECTA 2019 - 11th International Conference on Evolutionary Computation Theory and Applications

224

Figure 4: Hybrid CS-HS algorithm.

HS is widely used in medicine and engineering

design and in groundwater management problems

(Ayvaz, 2009). However, this method has a tendency

to population diversity decrease, which leads to

sticking at a local optimum.

Hybridization of the algorithms has the form of

sequential stages: each Cuckoo Search iteration is

followed by several Harmony Search iterations

applied to the decisions’ population. The right

adjustment of the Harmony Search properties has a

critical influence on the efficiency of the search. In

this work, we propose the particular settings of the

hybrid properties that were chosen after consideration

of the results of (Wang et al., 2016), and our own tests

on the multimodal analytic benchmark problems

(Jamil & Yang, 2013). First, HMCR was set to 1.0 to

make Harmony Search only to combine new solutions

by adjusting or remembering without improvisation.

This may quickly degrade the population, thus we set

the number of HS iterations not more than the size of

the population. This decision is empirical

experienced. HS and CS population sizes were set

equal. The values of all the properties are listed in

table 1.

Alongside the well-known HS/CS hybrid, two

improvements for it were considered (separately as

well as all at once). The first of them is also a well-

known modification by the author of the basic CS

Hybrid Cuckoo Search and Harmony Search Algorithm and Its Modiﬁcations for the Calibration of Groundwater Flow Models

225

method where the constant step size α is substituted

by chaotic varying step (Chaotic Cuckoo Search,

CCS) (Wang et al., 2016).

This significantly enhances the performance of

the CS by chaotically changing the decreasing

multiplier using chaotic sampling by a chaos map –

section, that was modified can be seen on, the upper

red box – where standard α is changed on chaotic α,

that is generated on each iteration. This hybrid is

referenced in the next section as CCS/HS.

More information on comparison of the

performance of the listed methods compared to the

classic heuristic algorithms like PSO and GA can be

found in (Wang et al., 2016), (Wang et al., 2016),

(Yang & Deb, 2009).

The modification to the HS component is

proposed by the authors of this work and its core idea

is to replace uniform random selection from the

memory with the probability HMCR by the selection

using modified normal sampling distribution. The

sampling steps are: 1) sample a normally distributed

value with 0 mean and HMS/3 (experience-based

value) standard deviation; 2) take absolute value and

discard non-integer part; 3) if this value is more than

the size of population – sample a value with uniform

distribution in [0, HMS-1] interval. The result value

will be the donor solution order number in population.

It allows us to use parameters of best solutions more

often while combining a new solution to evaluate

(population is sorted before constructing a new

solution). The hybrid method with modified HS is

referenced as CS/MHS. Modification of Harmony

Search is marked by the lower red box on.

And the last hybrid under consideration is the

combination of both chaotic CS and modified HS:

CCS/MHS.

Integral flow chart of described above hybrids is

presented on figure 4.

4 RESULTS

The basic hybrid method (CS/HS) and three proposed

variants were applied to the calibration of the two

groundwater flow models described above. The series

of 20 optimizations with random initial populations

were performed for each method and each model. The

target function for all cases was the sum of squared

residuals. The stop criterion for every optimization

was: 100 model runs without improvement of the best

found target value.

The model parameters, their ranges of variation

and the obtained optima results are given in the tables

2 and 3. Geological materials names are presented in

the tables 2 and 3 with their pseudonyms from the

figures 1 and 2. Models with optimized parameters

provided the agreement between the experimental

data and the simulation outputs at the observation

points with 5% precision. These metrics of the results

were quite similar for all four compared hybrids. The

obtained values allow making several conclusions.

First, the hydraulic conductivity coefficients for

element of geological structural model – parts of the

same lithology substantially vary with different

fracture density. It means that subselection of the

geological structural elements with diverse fracture

density may be crucial for the quality of simulations.

The second finding from the calibration is that

slightly fractured dykes are likely to act as natural

barriers.

The visualizations of the comparison of the

alternative algorithms via convergence plots

(Beiranvand, Hare and Lucet, 2017) are presented on

figures 5 and 6. The number of model runs is plotted

on the horizontal axis and the target function values

are on the vertical axis. The lower bound of each

filled contour is the convergence curve for the fastest

optimization among each method’s trials, the upper

bound – for the slowest one, and the bright-colored

line is the mean. According to the plots, all algorithms

appeared to be quite effective for model 1, their mean

convergence lines are very similar. Convergence for

model 1 is reached by 400-500 model runs. Model 2

required far more (about 10-20 thousand) model runs

to converge. And that is hardly surprising because of

more detailed parameterization of this model (7

parameters versus 3). And in this case, the

effectiveness of the considered algorithms

significantly differs. The CCS/HS hybrid converges

faster than others.

At first glance, the hybrids with modified

Harmony Search (CS/MHS and CCS/MHS) show

worse performance than CS/HS and CCS/HS on the

average.

On the series of convergence plots for these

algorithms (figures 7 and 8), one can see that

approximately half of the optimizations converges to

the suboptimal decision due to the tendency of the

modification to "condense" the decision set.

This ability to stably find not only the best

decision but also a set of suboptimal ones should be

considered further. In the context of the groundwater

flow and transport modelling the inherent uncertainty

of the available data is commonly rather significant.

And because of that suboptimal parametrization of

the model under consideration could be helpful for

the uncertainty assessment.

ECTA 2019 - 11th International Conference on Evolutionary Computation Theory and Applications

226

Table 2: Ranges of variation for model 1 parameters and found optimal values.

Material

Hydraulic conductivity coefficients [m/day]

Variation ranges

Optimal value

a1.Gneiss

(1.0∙10

-5

; 1,0∙10

-3

)

7.43∙10

-4

a2.Dikes

(1.0∙10

-5

; 1,0∙10

-3

)

6.93∙10

-4

a3.Quaternary sediments

(1.0∙10

-2

; 5,0∙10

-1

)

1.16∙10

-2

Table 3: Ranges of variation for model 2 parameters and found optimal values.

Material

Hydraulic conductivity coefficients [m/day]

Variation ranges

Optimal value

b1.Quaternary sediments

(1.0∙10

-3

; 5.0∙10

-1

)

4.998∙10

-1

b2.Weathering crust

(1.0∙10

-3

; 5.0∙10

-1

)

1.00∙10

-3

b3.Fissured dikes

(1.0∙10

-6

; 3.5∙10

–2

)

5.21∙10

-5

b4.Gneiss

(1.0∙10

-6

; 1.0∙10

-2

)

9.06∙10

-4

b5.Fissured gneiss

(5.0∙10

-3

; 1.0∙10

-1

)

5.001∙10

-3

b6.Сrushing zone

(1.0∙10

-2

; 1.0)

1.27∙10

-2

b7.Dikes

(1.0∙10

-6

; 1.0∙10

-3

)

6.09∙10

-2

Figure 5: Convergence plots for compared algorithms for

optimization, model 1.

Figure 6: Convergence plots for compared algorithms for

optimization, model 2.

5 CONCLUSIONS

In this work, the hybrid CS/HS optimization

algorithm and its modifications were successfully

Figure 7: Separate convergence plots for CS/MHS, model

Figure 8: Separate convergence plots for CCS/MHS, model

applied to the calibration of two alternative

groundwater flow models on the basis of geological

survey data of prospective site of the radioactive

waste deep geological disposal. The simulation

results obtained with the optimized parameters

Hybrid Cuckoo Search and Harmony Search Algorithm and Its Modiﬁcations for the Calibration of Groundwater Flow Models

227

appeared to be in good agreement with the

experimental data and could be used for radionuclides

transport simulations that are required as part of long-

term safety assessment for this sort of projects.

The combination of the proposed improvements

of the basic CS/HS algorithm was found to be

reasonable for this case of the optimization problem.

Hybrid methods with HMS component should be

developed further. And the CCS/HS variant appeared

to be the most efficient and stable among the others.

These qualities are highly valued for long-term safety

assessment purposes because the model calibration is

one of the key instruments (along with additional site

investigations) for the uncertainty treatment, and

accurate models are usually highly computationally

expensive. Hence, the proposed method could

become a noticeable contribution to the uncertainty

management framework within safety assessment

computational tools.

REFERENCES

Maier, H. R., Kapelan, Z., Kasprzyk, J., Kollat, J., Matott,

L. S., Cunha, M. C., ... & Ostfeld, A. (2014).

Evolutionary algorithms and other metaheuristics in

water resources: Current status, research challenges and

future directions. Environmental Modelling &

Software, 62, 271-299.

Haddad, O. B., Tabari, M. M. R., Fallah-Mehdipour, E., &

Mariño, M. A. (2013). Groundwater model calibration

by meta-heuristic algorithms. Water resources

management, 27(7), 2515-2529.

Fister, I., & Fister Jr, I. (Eds.). (2015). Adaptation and

hybridization in computational intelligence (Vol. 18).

Springer.

Hamid, O. H., & Braun, J. (2017, November).

Reinforcement Learning and Attractor Neural Network

Models of Associative Learning. In International Joint

Conference on Computational Intelligence (pp. 327-

349). Springer, Cham.

Wang, G. G., Deb, S., Gandomi, A. H., Zhang, Z., & Alavi,

A. H. (2016). Chaotic cuckoo search. Soft Computing,

20(9), 3349-3362.

Baldwin, T., Chapman, N., & Neall, F. (2008). Geological

disposal options for high-level waste and spent fuel.

Report for the UK Nuclear Decommissioning Authority.

Kapyrin, I., Konshin, I., Kramarenko, V., & Grigoriev, F.

(2018, September). Modeling Groundwater Flow in

Unconfined Conditions of Variable Density Solutions

in Dual-Porosity Media Using the GeRa Code. In

Russian Supercomputing Days (pp. 266-278). Springer,

Cham.

Jobmann, M. (2016). Site-specific evaluation of safety

issues for high-level waste disposal in crystalline rocks.

Final report (No. TEC--28-2015-AB). DBE Technology

GmbH.

Hill, M. C. (2004). Methods and guidelines for effective

model calibration. In Joint Conference on Water

Resource Engineering and Water Resources Planning

and Management 2000 (Vol. 104).

Wang, G. G., Gandomi, A. H., Zhao, X., & Chu, H. C. E.

(2016). Hybridizing harmony search algorithm with

cuckoo search for global numerical optimization. Soft

Computing, 20(1), 273-285.

Gandomi, A. H., Yang, X. S., & Alavi, A. H. (2013).

Cuckoo search algorithm: a metaheuristic approach to

solve structural optimization problems. Engineering

with computers, 29(1), 17-35.

Cole, B. J. (1995). Fractal time in animal behaviour: the

movement activity of Drosophila. Animal Behaviour,

50(5), 1317-1324.

Rhee, I., Shin, M., Hong, S., Lee, K., Kim, S. J., & Chong,

S. (2011). On the levy-walk nature of human mobility.

IEEE/ACM transactions on networking (TON), 19(3),

630-643.

Gupta, D., Das, B., & Panchal, V. K. (2013). Applying case

based reasoning in cuckoo search for the expedition of

groundwater exploration. In Proceedings of Seventh

International Conference on Bio-Inspired Computing:

Theories and Applications (BIC-TA 2012) (pp. 341-

353). Springer, India.

Ayvaz, M. T. (2009). Application of harmony search

algorithm to the solution of groundwater management

models. Advances in Water Resources, 32(6), 916-924.

Jamil, M., & Yang, X. S. (2013). A literature survey of

benchmark functions for global optimization problems.

arXiv preprint arXiv:1308.4008.

Laverov, N. P., Yudintsev, S. V., Kochkin, B. T., &

Malkovsky, V. I. (2016). The Russian strategy of using

crystalline rock as a repository for nuclear waste.

Elements, 12(4), 253-256.

Yang, X. S. (2010). Nature-inspired metaheuristic

algorithms. Luniver press.

Yang, X. S., & Deb, S. (2014). Cuckoo search: recent

advances and applications. Neural Computing and

Applications, 24(1), 169-174.

Beiranvand, V., Hare, W., & Lucet, Y. (2017). Best

practices for comparing optimization algorithms.

Optimization and Engineering, 18(4), 815-848.

Yang, X. S., & Deb, S. (2009, December). Cuckoo search

via Lévy flights. In 2009 World Congress on Nature &

Biologically Inspired Computing (NaBIC) (pp. 210-

214). IEEE.

ECTA 2019 - 11th International Conference on Evolutionary Computation Theory and Applications

228