Maximizing Sustainability in Reservoir Operation under Climate Change Using a Novel Adaptive Accelerated Gravitational Search Algorithm

Abstract

Holding a lasting balance between the water resources and water demands has become a challenging task for water resources managers, especially in recent years with the looming global warming crisis and its resulting climatic change effects. This paper focuses on modeling the optimized operation of the Zayandehrud Reservoir, located in west-central Iran, under two fifth-generation climate change scenarios called RCP4.5 and RCP8.5. A novel variant of the gravitational search algorithm (GSA), named the adaptive accelerated GSA (AAGSA) is proposed and adopted as the optimizer of the reservoir operation in this paper. The major advancement of the AAGSA against the original GSA is its high exploration capability, allowing the proposal to effectively tackle a variety of difficulties any complex optimization problem can face. The goal of the optimization process is the maximization of the sustainability of supplying the downstream water demands by the reservoir. The optimal results obtained by the original GSA and the proposed AAGSA algorithms suggest that the AAGSA can achieve much more accurate results with much less computational runtime, such that the proposed AAGSA is able to achieve the reservoir operation sustainability index of 98.53% and 99.46%, under RCP4.5 and RCP8.5 scenarios, respectively. These figures are higher than those obtained by the original GSA by 23.5% and 16% under RCP4.5 and RCP8.5, respectively, while the runtime of the proposal is reduced by over 80% in both scenarios, as compared to the GSA, suggesting the high competence of the proposed AAGSA to solve such a high-dimensional and complex real-world engineering problem.

1. Introduction

In recent years, climate change has turned into a critical phenomenon as it takes considerable effects on decreased available water resources. Increased population, urban construction, land erosion, land-use change, and mismanagement of the water resources have all made the climate change effects doubly important and more necessary to mitigate [1,2]. The variations in the time, intensity, and type of precipitation, the increased frequency of the extreme phenomena such as floods and droughts, the decreased snowfall, and the intensified ice melting, are among the major and tangible consequences of climate change [3,4,5]. These events may take negative and irreversible effects on human health and the environment if no planning is suggested to adapt the water management policies to climate change [6]. The research shows that climate change could pose significant changes in the inflow and the release of the surface water reservoirs. Since the reservoirs have mainly been designed and operated to achieve a wide range of purposes such as supplying the domestic, industrial, and agricultural water demands, generating the hydro-power energy, controlling the floods, and navigation and tourism purposes [7], the effects of climate change on the reservoirs could pose increased conflict and competition between the stakeholders and water consumers, deepening the water crisis and thus, highly urging the managers to revisit the reservoir operation policies even more than before [8,9]. As a result of the increased water demands, conflicts between users of any water resources will arise and this is the more valuable water users’ demand which would have the decision-makers consider reallocation policies to avoid over-exploitation of water resources [10].

In many of the research studies recently conducted on reservoir operation under climate change, numerous optimization methods along with several climate change scenarios have been utilized to better project the climatic components and the resulting inflow to the reservoirs [11,12,13,14]. All these methods and scenarios are based on the general circulation models (GCMs) to simulate the global climate and make long-term projections using the Intergovernmental Panel on Climate Change (IPCC) emission scenarios as the boundary conditions of these models. The GCM models also need the downscaling methods to project the climate variables at the regional and sub-regional scales [15]. In the reservoir operation problems, the precipitation and temperature are first projected by the GCM models. Thereafter, they are downscaled and fed to a rainfall-runoff processing model to yield the inflow to the reservoirs. For example, Zamani et al. [16] utilized LARS-WG to downscale the outputs of 14 GCM models under the A2 and B1 scenarios and entered the results into the IHACRES rainfall-runoff model to obtain the inflow to a reservoir in Iran. The purpose was to optimize the performance of the reservoir. In this research, projecting the climate variables utilizing several GCM models and weighting these models using the sensitivity analysis is found to be more accurate than projecting the climate variables using only one GCM model. Besides the use of accurate projection models, adopting the appropriate optimization method is also necessary to optimize the reservoir operation under climate change conditions [17]. Various mathematical optimization methods have been applied to reservoir operation optimization problems. Dobson et al. [18] reviewed the literature of the works carried out in the field of reservoir operation optimization and presented a novel classification system focusing on the different arguments of the optimization problems in this specific field. Recently, the meta-heuristic optimization algorithms open a new horizon to effectively and reliably solve the intrinsically complex, high-dimensional, and hard-to-handle reservoir operation optimization problems [19,20]. In most recent research studies, various optimization algorithms are developed to face complex problems. Rezaei et al. [21] proposed the VAGWO algorithm in which a velocity term is added to the position-updating procedure of the original GWO algorithm to improve the exploration capability of GWO. The results proved it’s the proposed algorithm can yield higher-quality results with a plausible computational cost compared to GWO in solving a variety of high-dimensional and complex benchmark functions and several engineering problems. To solve a reservoir operation problem, many input parameters should be inserted into the optimization model such as precipitation, inflow, evaporation from the reservoir surface, and the downstream demands. The output of such optimization models could generally be the optimal water releases in different time steps. Jothiprakash and Shanthi [22] utilized the genetic algorithm (GA) to minimize the square water shortages while holding the reservoir storage at a desirable level. Cheng et al. [23] introduced the chaotic genetic algorithm (CGA) in which the chaotic operators replace the ordinary GA operators. The results suggested that the CGA could properly avoid local optima and increase the hydro-power energy generated by the reservoir and decrease the convergence rate of the algorithm, as compared to the dynamic programming (DP) and the original GA. Rezaei and Safavi proposed DFPSO algorithm by integrating the abilities for evaluation of fitness and diversity of the particles in the PSO algorithm to avoid unbalanced solutions’ movements in determining their global guides. After being verified by several benchmark functions, the proposal was applied to a conjunctive surface–ground water use problem to minimize the shortages in meeting irrigation water demands under different climatic conditions while holding the groundwater resources sustainable to use. The results showed its high efficacy and efficiency in solving such a complex engineering problem [24]. Zhang et al. [25] proposed an improved particle swarm optimization (PSO) algorithm, named IAPSO, which can dynamically adjust the cognitive and social scaling parameters to avoid premature convergence. The proposed algorithm was then applied to a hydro-power reservoir operation optimization problem and demonstrated to be an efficient algorithm to boost the hydro-power generation with a high convergence speed. Asvini and Amudha [26] applied the plant propagation algorithm (PPA) to optimal reservoir operation in a case study in India and the results suggested the high performance of the algorithm to meet a large portion of the water demands even in the dry months of the planning period. Bozorg-Haddad et al. [27] examined the capability of the GSA to solve the optimization problems via the implementation of this algorithm on three benchmark functions, a hydropower reservoir, and a four-reservoir optimization problem. The results suggested the ability of the GSA to faster and also more accurately converge to the optimal solution as compared to the GA. Turgut et al. [28] introduced the master–slave and crow search algorithms and applied them to the long-term operation of a high dam reservoir. The results revealed the superiority of the master–slave algorithm to solve this problem mainly due to its high ability to realize the uncertain and non-linear nature of the reservoir systems.

Recently, hybrid models and algorithms have been developed and used in various research studies to produce more reliable outcomes. Alizamir et al. [29] investigated two hybrid models ANFIS-PSO and ANFIS-GA to calculate reference evapotranspiration (ET₀) in Antalya and Isparta stations in Turkey. The results showed that both hybrid models are able of generating more accurate monthly ET₀ than those resulting from ANN, ANFIS, and CART. Kadkhodazadeh et al. [30] developed the LSSVM-GBO algorithm to assess water quality parameters EC and TDS in Ahvaz, Armand, and Gotvand stations in the Karun River basin after validating their performance via applying them to three benchmark datasets. The results revealed the LSSVM-GBO to be superior to the ANN and ANFIS.

Several different indices are already proposed in the literature to benchmark the performance of the water resources systems including: (1) drought severity index [31]; (2) environmental sustainability index [32]; (3) integrated drought index [33]; (4) water system performance index (WSPI) [34]; and (5) sustainability index [35]. The sustainability index is a combination of three well-known and widely used criteria including reliability, resilience, and vulnerability (RRVs) [36]. Sandoval-Solis et al. [37] modified the sustainability index by improving the structure and scale of this index to be well adjusted to any water use sector and any water basin throughout the world. There are numerous methods proposed in the literature to more reliably and more accurately calculate the RRVs index [38,39,40]. Safavi and Golmohammadi [41] proposed a new approach utilizing fuzzy logic to calculate the sustainability index of the Zayandehrud Reservoir operation in west-central Iran. They compared their proposed approach with those calculating the sustainability index through the classical equations. They concluded that the fuzzy index could take into account the uncertainties in meeting a portion of water in the temporal intervals and allow the experts and the stakeholders to have more tangible data ahead and enable them to more reliably judge on the level of sustainability held in the water resources systems operation. Several other research studies have been carried out on the role of fuzzy logic to improve the notion of the sustainability of the water resources systems in the literature [42,43,44].

The Zayandehrud River basin is located in the arid and semi-arid climatic conditions in west-central Iran. Since a streak of drought events has occurred in this river basin, especially in recent years, the optimal management of the only large dam operating in this basin is of high importance [45,46,47]. Accordingly, this reservoir is investigated as a case study of this paper to be optimally managed when operating to meet the water demands in a near-future climate-change-affected period.

Figure 1 illustrates a brief flowchart of the research carried out in this paper. As can be seen in the figure below, the evaporation from the surface of the reservoir is first estimated using the Torrent White method, based on the daily temperature data processed and downscaled by the LARS-WG model. These data along with the precipitation data are already inserted into the HEC-HMS rainfall-runoff model under the RCP4.5 and RCP8.5 climate change scenarios to simulate the inflow to the reservoir. All of these preliminary stages are already accomplished by Sangestani [48], Shamaeian [49], and Alkantar [50]. The novelty of this paper is in proposing a new variant of the gravitational search algorithm (GSA), named adaptive accelerated GSA (AAGSA) to maximize the sustainability in a reservoir operation problem. The main objective of the optimization process is to minimize the water shortages to meet the total downstream demands while preserving the reservoir storage sustainability to enable the stakeholders to consistently operate the reservoir in the years following the planning period. The performance of the AAGSA is then compared to the original GSA to benchmark the best-performing algorithm to yield the maximum sustainability index as a criterion explaining the degree of satisfaction for the demands while controlling the water stored in the reservoir.

The remainder of this paper is organized as follows. Section 2 explains how to prepare climate data under the climate change scenario and estimate the inflow to a reservoir as the main characteristic of a reservoir system affected by climate change. In the following, this Section describes the methodology of this paper, reviewing the gravitational search algorithm (GSA) and proposing the AAGSA algorithm developed to improve the key capabilities of the GSA. Section 3 introduces the study area and presents the optimization model framework to be solved in this paper. Section 4 presents the results and discusses them. Finally, Section 5 concludes the paper.

2. Materials and Methods

2.1. Projection of the Future Climate Variables

In the last few years of the 20th century, the evidence of the effect of global warming on the earth’s climate system, contributed to the founding of the Intergovernmental Panel on Climate Change (IPCC). The mission of this panel is to scientifically evaluate the climatic changes throughout the world, the resulting consequences of these changes, and propose effective and practical solutions to decrease the likely dangers of climate change that threaten the lives of humans and the other living bodies in all over the world. According to the recent reports of the IPCC, the mean temperature of the world in 2017 has reached 1 °C more than that recorded in the Industrial Revolution era, and it may even approach 1.5 °C above than that in the Industrial Revolution era until 2040, imposing irreversible damages to humankind and the environment. Following the negative consequences of climate change, the IPCC mandates countries to decrease the emission of greenhouse gases by the extent adjusted to the share of each country to emit these gases to the atmosphere [6]. The IPCC proposes several climate change scenarios to approximately estimate the effect of climate change in the long-term future period to enable the policymakers to better adjust the development schemes to the climatic conditions ahead. The representative concentration pathway (RCP) scenarios are first proposed in 2009 and prescribed by the IPCC to climatically model the world in 2014. These scenarios fall into four main categories named RCP2.6, RCP4.5, RCP6, and RCP8.5, based on the level of the radiation induction resulting from the greenhouse gases emitted to the atmosphere until 2100. As a result, the RCP2.6 is taken as the most optimistic scenario and the RCP8.5 scenario is rated as the most pessimistic one.

To project the temporal events occurring in the atmosphere, the general circulation models (GCMs) are usually used. These models can solve the equations governing the atmosphere to simulate the earth’s climate within spatial intervals. The dimensionality of the spatial simulations performed by these models is 200 km by 200 km. The GCMs use the boundary conditions assumed in the IPCC scenarios to make long-term projections on the climate variables. As a result, the GCMs are unable to directly project the climate variables in the pointwise or regional scales, but their projected data must be downscaled to fit any desired scale [51]. There are two methods to carry out downscaling: (1) dynamical methods and (2) statistical methods. In this paper, the long Ashton Research Station Weather Generator (LARS-WG) model is used as a statistical model to estimate the climatic outputs with an accuracy of 0.5 degrees. In the statistical downscaling, the univariate or multi-variate regression, artificial neural networks (ANNs), or some other statistical approximation techniques are utilized to hold a relation between the real behavior of the climate station and the output of the general circulation models (GCMs). The statistical methods are more interested among the researchers, as they are more computationally efficient and also easier to implement compared to the other downscaling methods. The LARS-WG was first developed in 1991 in an agricultural risky evaluation project conducted by the Hungarian scientists’ academy. This model can generate a time series of the daily climatic data such as the precipitation, maximum and minimum temperature, and the solar radiations as the outputs of the model. The inputs of the LARS-WG include the characteristics of the climate stations such as their names and the coordinates, and also a list of the observed daily climatic data. The data generation process conducted by this model is summarized in three main steps as follows:

• Calibration (site analysis): analyzing the observed data and saving the information reached in two separate folders.
• Validation (Qtest): comparing the data artificially generated by the model with the data observed and inserted to the LARS-WG based upon the substantial statistical differences between these two data sets.
• Generation: producing a time series of the spatially downscaled future temporal data with a daily scale [52].

2.2. Runoff Estimation

To estimate the runoff occurring in the future planning period, the HEC-HMS rainfall-runoff model is utilized in this paper. The HEC-HMS is one of the hydrological modeling software developed to simulate the rainfall-runoff process in the dendritic basins systems. This model is able to divide the hydrological cycle into several controllable pieces and discriminating the boundaries in the river basin to accurately build up the rainfall-runoff processing model. Each mass or energy in this model can be represented by a mathematical model. The main sections of the HEC-HMS include the river basin’s physical characteristics, climatic characteristics, hydrological simulations, parameter tuning, calibration stage, and finally validation stage. The latter stage is to evaluate the degree of adaptation between the computed data and the observed data regarding a number of the performance criteria existing in the model [53].

In this paper, the HEC-HMS rainfall-runoff model was fed by the inputs including the precipitation, and maximum and minimum daily temperature data obtained for the six selected stations in the upstream basin of the Zayandehrud Reservoir. These inputs are the outputs of the LARS-WG model that were achieved and reported in [48,49], after the model was built up, calibrated, and validated in [50]. The HEC-HMS was then used to generate the inflow to the reservoir for the 2020–2033 future planning period under the RCP4.5 and RCP8.5 scenarios. Since there are no significant and tangible differences between the outputs of the simulation models under either of the RCP scenarios in the near-future period, the RCP4.5 scenario was adopted as a medium-standing scenario, representing all of the other scenarios to be used in the reservoir operation optimization process under climate change conditions. The reservoir inputs under the RCP8.5 scenario were also chosen to generate input data for the HEC-HMS simulation model. The simulated inflow data provided were then imported to the optimization models. Thereafter, the outcomes of the optimization models including the sustainability index were evaluated to be compared with those offered by the models run by the data produced by the RCP4.5. The inter-basin water transferred through the tunnels in the future planning period was also simulated by the WEAP model in [48]. The water volume transferred was added to the natural inflow to the reservoir to form the total inflow estimated under the climate change scenarios for the planning period in this study.

2.3. Gravitational Search Algorithm

The gravitational search algorithm (GSA), first proposed by Rashedi et al. [54], is one of the well-known meta-heuristic optimization algorithms modeling the gravitational laws and the Newtonian motion in the framework of an artificial system operating in a discrete-time domain to search for the global optimum of the single-objective optimization problems. The GSA could be categorized in the metaphor-based, memory-free, and population-based algorithms, with a multi-neighborhood structure for the search agents that are aimed at optimizing a static objective function. In the GSA system search, each search agent has a mass calculated to reflect its fitness in the search space. Each search agent attracts another agent with an acceleration component directly proportional to its mass and inversely proportional to the distance between these two masses. In this algorithm, the search agents with the best fitness are called the elite agents. The number of these elite agents is high at the early stages of the optimization process and is gradually decreased when approaching the final stages of this process.

Generally, the more the mass of an elite search agent, the more the tendency of a search agent to move toward that elite agent provided that the distance between the two agents is short enough. When an agent faces an elite one having a long distance to it, the agent deaccelerates its motion toward the elite agent provided that the mass of the elite agent is low enough to justify this deacceleration made to the movement. The main reason for this behavior of the search agents in that the search space may be hidden in the exploration and exploitation mechanisms considered to be held in the GSA system search. When there is a long distance between an elite agent and an ordinary one, there are a lot of points undetected in the search space and the neighborhood of the ordinary agent, urging that agent to retard its movement through the elite agent to become able to find out any further good positions in this course. Nevertheless, if the mass of the elite agent is large enough, the ordinary agents could get more accelerated when moving toward the elite agent. In this way, both the exploration and the exploitation processes may be better accomplished during the optimization process. Furthermore, there might be a smooth and effective exploration–exploitation balance in the search process through gradually decreasing the number of the elite agents by lapse of the optimization stages. The stages of the optimization could be modeled as the discrete iterations by lapse of which the heaviest agent (the fittest point) in the search space may be found before reaching the pre-defined stopping criterion of the optimization process.

Suppose for an N-agent system. The position of the ith agent at the dth dimension may be expressed as follows:

$$X_i=\left(x_i^1{,\dots ,x}_i^d{,\dots ,x}_i^n\right) \mathrm{for} i=1, 2,\dots ,N$$

(1)

The gravitational force applied to the ith agent by the jth one at the tth iteration is calculated through Equation (2).

$$F_{\mathit{ij}}^d\left(t\right)=G\left(t\right)\frac{M_{\mathit{pi}}{\left(t\right)\times M}_{\mathit{aj}}\left(t\right)}{R_{\mathit{ij}}\left(t\right)+\epsilon }{(x}_j^d{\left(t\right)-x}_i^d\left(t\right))$$

(2)

$$R_{\mathit{ij}}\left(t\right){=||X}_i{\left(t\right),X}_j{\left(t\right)||}_2$$

(3)

where $M_{\mathit{aj}}$, and $M_{\mathit{pi}}$ are the active and passive masses of the jth and ith agents, respectively; $G\left(t\right)$ is the gravitational constant at the tth iteration; and $R_{\mathit{ij}}\left(t\right)$ is the Euclidean distance between the ith and jth agents at the tth iteration. $\epsilon$ is a small positive constant applied to hinder singularity. $G\left(t\right)$ can be calculated as follows:

$${G\left(t\right)=G}_0{e}^{-\frac{\alpha t}{T}}$$

(4)

where $G_0$ is the initial constant value which should be precisely set for $G\left(t\right)$, and $\alpha$ is a constant which should be carefully set to significantly reduce the $G\left(t\right)$ at the final iterations to allow the exploitation process to be started. T is the total number of iterations representing the life span of the system. The resultant force applied to the ith agent by the elite agents can be calculated as shown in Equation (5).

$$F_i^d\left(t\right)={\displaystyle \sum }_{j\in \mathit{Kbest} , j\ne i}{\mathit{rand}}_j^d{F}_{\mathit{ij}}^d\left(t\right)$$

(5)

$$\mathit{Kbest}=\left({\mathit{final}}_{\mathit{per}}+\left(1-\frac{t}{T}\right)\times \left({100-\mathit{final}}_{\mathit{per}}\right)\right)\times \frac{N}{100}{, \mathit{final}}_{\mathit{per}}=2$$

(6)

where $Kbest$ is the dynamic number of the elite agents adjusted with regard to the number of iterations as formulated in Equation (6). $ran{d}_j^d$ is a uniformly distributed random number that should be generated in [0, 1] and assigned to the dth dimension of the jth agent. These random numbers are mainly generated in the GSA to impart a stochastic nature to any movement made by the search agents in this algorithm. Finally, the acceleration the ith search agent takes in the dth dimension at the tth iteration denoted by $a_i^d\left(t\right)$ can be calculated as follows:

$$a_i^d\left(t\right)=\frac{F_i^d\left(t\right)}{M_{\mathit{ii}}\left(t\right)}$$

(7)

where, $M_{\mathit{ii}}\left(t\right)$ is the inertial mass of the ith agent which is assumed to be equal to the active and passive masses of the corresponding agent in the original version of the GSA. The updating equation for the velocity and the position of each agent in the search system of the GSA is as follows:

$$v_i^d\left(t+1\right){=\mathit{rand}}_i{\times v}_i^d{\left(t\right)+a}_i^d\left(t\right)$$

(8)

$$x_i^d\left(t+1\right){=x}_i^d\left(t\right){+v}_i^d(t+1)$$

(9)

where, $v_i^d\left(t+1\right)$ and $v_i^d\left(t\right)$ denote the updated and the current velocity of the ith agent; $x_i^d\left(t+1\right)$ and $x_i^d\left(t\right)$ are the updated and the current position of the ith agent; $a_i^d\left(t\right)$ is the acceleration the ith agent takes at the current iteration; and ${\mathit{rand}}_i$ is a random number in [0, 1] generated and assigned to the current velocity of an agent. Since finding a better agent with the heavier mass in the GSA memory-less search system is not guaranteed at any iteration, applying a decreasing inertial weight in this algorithm is avoided and instead, these random numbers are employed to gradually make the agents static in the search space.

The masses of the agents at each iteration are calculated according to their fitness as follows:

$$m_i\left(t\right)=\frac{{\mathit{fit}}_i\left(t\right)-\mathit{worst}\left(t\right)}{\mathit{best}\left(t\right)-\mathit{worst}\left(t\right)}$$

(10)

$$M_i\left(t\right)=\frac{m_i\left(t\right)}{{{\displaystyle \sum }}_{j=1}^N{m}_j\left(t\right)}$$

(11)

$$\mathit{best}\left(t\right)=\underset{j\in \{1,\dots ,N\}}{\mathit{min}}{\mathit{fit}}_j\left(t\right)$$

(12)

$$\mathit{worst}\left(t\right)=\underset{j\in \{1,\dots ,N\}}{\mathit{max}}{\mathit{fit}}_j\left(t\right)$$

(13)

where ${\mathit{fit}}_i\left(t\right)$ is the fitness function of the optimization problem. Furthermore, Equations (12) and (13) are set for the minimization problems. For the maximization problems, the min and max in Equations (12) and (13) must be replaced with max and min, respectively.

Equation (7) can be rewritten according to Equation (5) and after some simplifications, as follows:

$$a_i^d\left(t\right)=G\left(t\right)\times {\displaystyle \sum }_{j\in \mathit{Kbest} , j\ne i}{\mathit{rand}}_j^d\times \frac{M_j\left(t\right)}{R_{\mathit{ij}}\left(t\right)+\epsilon }\times \left(x_j^d\left(t\right){-x}_i^d\left(t\right)\right)$$

(14)

The steps of the GSA are as follows [17]:

• Defining the bounds of the search space.
• Randomly generate the positions of the search agents.
• Fitness evaluation of the search agents.
• Updating the $G\left(t\right), \mathit{best}\left(t\right), \mathit{worst}\left(t\right), M_i\left(t\right), \mathrm{for} i=1, 2, \dots , N$.
• Calculating the resultant force applied to every agent in all dimensions.
• Calculating the acceleration of all agents.
• Updating the velocity and position of the agents.
• Repeating steps 3 to 7 until reaching the stopping criterion.
• Stop.

2.4. Proposed Algorithm: Adaptive Accelerated Gravitational Search Algorithm

The main weakness of the GSA may be found in its procedure to calculate the acceleration of the search agents as the major novelty discriminating this algorithm from the others. Accordingly, an improved mechanism is proposed in this paper to better define an acceleration component for each search agent to impede the premature convergence for the GSA. This premature convergence mainly occurs when the original GSA generally holds the acceleration components too small and thus, slows down the global search procedure and disrupts the exploration phase of the optimization process. This drawback of the GSA may be exacerbated when the optimization problem the GSA attempts to solve is of high complexity and/or high dimensionality. To alleviate the drawbacks mentioned above and to enhance the exploration capability of the GSA, a novel variant of the GSA, named adaptive accelerated gravitational search algorithm (AAGSA), is proposed in this paper. The main improvements the AAGSA makes to the original GSA can be summarized in two major points:

• In the AAGSA, the distances between each couple of the search agents are also normalized and transformed to the values lying in [0, 1], as the masses remain normalized. This procedure can, on the one hand, better show the differences among several distances an elite search agent may have with the ordinary agents attempting to move toward it, and on the other hand, can impede these distances from rapidly growing, especially in the high-dimensional problems. In this way, the exploration capability could be boosted, as the value of the acceleration term of the updating procedure can be significantly increased to better solve the optimization problems, especially the high-dimensional and complex problems. The secondary effect of the normalization of the distances in the AAGSA is rectifying the dimensionality of the acceleration term. As can be seen in Equation (8), the acceleration term would be added to the velocity term. As a result, the dimensionality of the acceleration is necessary to be of the length dimensionality, whereas in the original GSA, the acceleration term is dimensionless. While the normalization of the distances between the agents makes the acceleration be of length dimensionality. In this way, this technical drawback of the GSA is removed.
• In the original GSA, all the random numbers multiplied by the forces are generated in [0, 1], disregarding if the algorithm is in the early or the late iterations of the optimization process. While a more efficient exploration–exploitation balance can be held by the proposed AAGSA algorithm through emphasizing the masses of the agents at the early iterations less, while emphasizing the distances in these iterations more. By lapse of iterations, the influence of the masses is gradually increased and the influence of the distances is decreased. These modifications can be justified as the masses are assumed to be of more uncertainty at the early iterations and are deemed to be of less uncertainty at the later iterations. On the other hand, the distances are assumed to have more effect in guiding the agents at the early iterations while their effect gradually decreases and is equated to the effect of the masses when approaching the final iterations. These adjustments are mainly made to the AAGSA to help it diversify the search agents at the initial iterations and intensify the local search in achieving the high-fitness areas in the search space at the final iterations. These adjustments can be made to the AAGSA via modifying the ranges the random numbers are to be generated in, such that the lower bounds and the upper bounds of these ranges are dynamically growing from the upper neighborhood of 0 to the lower neighborhood of 1 in the general interval [0, 1]. In this procedure,${ c}_1\left(t\right) \le { \mathit{rand}}_{j }^d\le { c}_2\left(t\right)$, and ${\mathit{rand}}_j^d$ is the random number generated for the jth attracting agent at the dth dimension, while the parameters $c_1\left(t\right)$ and $c_2\left(t\right)$ can be calculated based on Equations (15)–(19). The resulting ranges are thus made something similar to [0, 0], [0, 0.1], [0.1, 0.2], [0.2, 0.3], …, [0.9, 1], and [1, 1] as the iterations go on.

  $$c_1\left(t\right)=1-\frac{1-\frac{t}{T}}{1-{\mu }_1\left(t\right)\times \frac{t}{T}}$$

  (15)

  $$c_2\left(t\right)=1-\frac{1-\frac{t}{T}}{1-{\mu }_2\left(t\right)\times \frac{t}{T}}$$

  (16)

  $${\mu }_1\left(t\right)=1-{\left(\frac{t}{T}\right)}^{\mathit{power}\left(t\right)}$$

  (17)

  $${\mu }_2\left(t\right)=1-\frac{t}{T}$$

  (18)

  $$\mathit{power}\left(t\right){=\mathit{power}}_{\mathit{max}}-{(\mathit{power}}_{\mathit{max}}{-\mathit{power}}_{\mathit{min}})\times \frac{t}{T}$$

  (19)

In these equations, $\mathit{power}\left(t\right)$ is a control parameter linearly decreasing from ${\mathit{power}}_{\mathit{max}}$ to ${\mathit{power}}_{\mathit{min}}$ as the iterations progress. ${\mathit{power}}_{\mathit{max}}$ is set to 4 and the ${\mathit{power}}_{\mathit{min}}$ is set to 1 in this paper. During the optimization process, $c_1\left(t\right)$ and $c_2\left(t\right)$ consistently increase to make both the upper bound and lower bound of the ranges of the random numbers go ahead and reach [1, 1] = $\left\{1\right\}$. The role of the parameter $\mathit{power}\left(t\right)$ can be found to make ${\mu }_1\left(t\right) \ge { \mu }_2\left(t\right)$ and so, $c_1\left(t\right) \le { c}_2\left(t\right)$. The parameters ${\mu }_1\left(t\right)$ and ${\mu }_2\left(t\right)$ have the mission of making the parameters $c_1\left(t\right)$ and $c_2\left(t\right)$ larger than what is normally expected, contributing the ranges of the random numbers to be consistently closer to [1, 1] = $\left\{1\right\}$ than what is normally expected. These mathematical properties based on which the random numbers are frequently generated in AAGSA can be very effective in strengthening the exploitation capability. This improvement made to the GSA along with enhancing the exploration capability of this algorithm which is already carried out through normalizing the distances could comprehensively strengthen the GSA in the framework of the proposed AAGSA. Accordingly, Equations (20)–(22) propose a new definition for the acceleration term as is set in the AAGSA.

$${\mathit{rand}}_j{=c}_1\left(t\right){+\mathit{rand}}_c\times \left(c_2\left(t\right){-c}_1\left(t\right)\right)$$

(20)

$$a_i^d\left(t\right)=G\left(t\right)\times {\displaystyle \sum }_{j\in \mathit{Kbest} , j\ne i}{c}_{\mathit{random}}{\times \mathit{rand}}_j\times \frac{M_j\left(t\right)}{\left(R_{\mathit{ij}}\left(t\right){/R}_{\mathit{max}}\right)+\epsilon }\times \left(x_j^d\left(t\right){-x}_i^d\left(t\right)\right)$$

(21)

$$R_{\mathit{max}}=\sqrt{{\displaystyle \sum }_d{{(x}_{\mathit{max}}^d{-x}_{\mathit{min}}^d)}^2}$$

(22)

where ${\mathit{rand}}_c$ is a random number generated in [0, 1] to make the ${\mathit{rand}}_j$ be produced in the range set $\left[c_1\left(t\right){, c}_2\left(t\right)\right]$. $c_{\mathit{random}}$ is a coefficient set to 2 to further help the agents not be trapped in local optima. $x_{\mathit{max}}^d$ and $x_{\mathit{min}}^d$ are the upper and lower numerical values of the search space at the dth dimension of this space, and $R_{\mathit{max}}$ is the maximum distance between any two points in the search space to be used to normalize the distances represented by $R_{\mathit{ij}}\left(t\right)$ in Equation (21). Finally, $a_i^d\left(t\right)$ is the acceleration term redefined in the proposed AAGSA and used in the updating procedure of this algorithm. It is noteworthy that since the diversification process of the AAGSA is set to be much better than that of the GSA, the random numbers generated to make the search process stochastic are decided to be varied by only the agent and are constant when the dimension varies. This is why this random number is represented by ${\mathit{rand}}_j$ rather than ${\mathit{rand}}_j^d$ in Equations (20) and (21). This modification can also significantly reduce the computational cost of the AAGSA compared to that of the GSA, especially when handling high-dimensional optimization problems.

The secondary modification applied to the GSA in the framework of the AAGSA is to redefine the bound constraint handling of this algorithm. While in the GSA an infeasible variable is randomly reinitialized between the upper and lower bounds based on a uniform probability distribution, a uniform mutation is imposed on these infeasible variables in AAGSA. In this mutation, the mutation coefficient is large at the early iterations and is dwindling by lapse of iterations to further keep the agents in the search space at the early stages of the optimization and to gradually decrease pushing the infeasible agents to be feasible when approaching the later stages. Equations (23) and (24) show the constraint handling mechanism in the original GSA.

$$x_i^d\left(t+1\right){=x}_{\mathit{min}}^d+\mathit{rand}\times \left({ x}_{\mathit{max}}^d{-x}_{\mathit{min}}^d\right){ ; \mathrm{if} x}_i^d\left(t+1\right){<x}_{\mathit{min}}^d$$

(23)

$$x_i^d\left(t+1\right){=x}_{\mathit{max}}^d-\mathit{rand}\times \left({ x}_{\mathit{max}}^d{-x}_{\mathit{min}}^d\right){ ; \mathrm{if} x}_i^d\left(t+1\right){>x}_{\mathit{max}}^d$$

(24)

Equations (25)–(27) illustrate the constraint handling procedure proposed in the AAGSA.

$$x_i^d\left(t+1\right){=x}_{\mathit{min}}^d+\mathit{rand}\times \mathit{mut}(t+1)\times \left({ x}_{\mathit{max}}^d{-x}_{\mathit{min}}^d\right){ ; \mathrm{if} x}_i^d\left(t+1\right){<x}_{\mathit{min}}^d$$

(25)

$$x_i^d\left(t+1\right){=x}_{\mathit{max}}^d--\mathit{rand}\times \mathit{mut}(t+1)\times \left({ x}_{\mathit{max}}^d{-x}_{\mathit{min}}^d\right){ ; \mathrm{if} x}_i^d\left(t+1\right){>x}_{\mathit{max}}^d$$

(26)

$$\mathit{mut}\left(t+1\right){=\mathit{mut}}_{\mathit{max}}-\left({\mathit{mut}}_{\mathit{max}}{-\mathit{mut}}_{\mathit{min}}\right)\times \frac{t}{T}{; \mathit{mut}}_{\mathit{max}}=0.9; {\mathit{mut}}_{\mathit{min}}=0.1$$

(27)

In the original GSA, there is no constraint on the velocities of the agents at any dimension. Hence, as the final modification applied to the AAGSA, for the velocity of the search agents in all dimensions an upper bound and a lower bound are also defined, exceeding which is assumed as the sign of infeasibility of the decision variables represented by the search agents’ positions in the search space.

3. Study area

3.1. Zayandehrud Dam Reservoir

The Zayandehrud Dam is constructed on the Zayandehrud River as the largest river running in the Isfahan Province at the central plateau of Iran (Figure 2). The Zayandehrud River basin is the most Zayandehrud Dam-consuming sector and is located between 50°24ʹ–53°24ʹ longitudes and 31°11ʹ–33°42ʹ latitudes. Covering an area of 26,917 km², this river basin is rated the most important basin in central Iran necessitating the Zayandehrud Dam as the main surface water storage system upstream of this basin to supply a large variety of domestic, agricultural, industrial, and environmental water demands of the basin [55]. The Zayandehrud Reservoir operates to supply the drinking and sanitary water for 5.2 million people residing in three provinces included in the Zayandehrud basin. In addition, this reservoir is in charge of supplying water demands of more than 200,000 ha of the agricultural areas in the basin. Furthermore, most of the largest industries of Iran are centralized and located in this basin such as the steel and cement production factories, power stations, oil refineries, and so on. A large portion of the water requirements of these industries is supplied by the Zayandehrud Reservoir. Moreover, since the Zayandehrud River, as one of the longest and most important Rivers in Iran, is running through the globally reputed and tourist city of Isfahan located at the center of the Zayandehrud Basin, this Basin is of high importance in terms of the tourism industry and the relevant industries that can grow and flourish due to the tourism industry. As a result, the perennial water stream in the Zayandehrud River is one of the other factors of high importance that can be supported by the proper and adequate water released from the Zayandehrud Reservoir. Moreover, the water resources of this basin are transferred to the adjacent basins to supply the water of two other large cities of Iran for drinking and sanitary purposes. In the past decade, in particular, the basin was suffering from a streak of severe droughts accompanied by increased population and thus, increased water demands in the domestic and agricultural sectors. All of these characteristics counted for the Zayandehrud basin have made an imbalance between the water resources, especially the surface water resources stored in the Zayandehrud Reservoir, and the water demand by all sectors. In recent years, global warming and its resulting climatic changes have significantly escalated the imbalance between these water resources and water demands and have strongly encouraged the water managers, planners, and decision-makers to adapt their plans to the newly emerging and specific conditions for the water resources and water demands in this important basin, as much as possible. As a result, the Zayandehrud basin has been made a very strategic region needing the consistent attention of the authorities when planning the water resources management schemes for the basin. The points expressed in this sub-section could suitably unravel the importance of the Zayandehrud Reservoir to be addressed as the case study of the reservoir operation in this paper [55]. The water inflowing to this reservoir is supplied by the direct precipitation, three rivers in the basin, and also three tunnels transferring the water from the adjacent basin of this basin [56,57]. The location of the Zayandehrud Dam and Reservoir in the Zayandehrud upstream basin is displayed in Figure 2.

Table 1. The detailed characteristics of the Zayandehrud Dam and Reservoir.

Characteristics	Statements
Dam type	Arch dam
Coordinates	32.73° N
	50.74° E
Total capacity	1470 MCM
Dead storage	150 MCM
Average annual inflow	About 1400 MCM
Dam lake area	54 km2

summarizes the detailed characteristics of the Zayandehrud Dam and Reservoir.

As mentioned in the reports, the Zayandehrud Reservoir must operate to supply 400 million cubic meters (MCM) annually as the domestic (drinking and sanitary) water demands, while the environmental water needs of the Gavkhouni wetland is estimated to be 176 MCM, which is never fully supplied by the Zayandehrud Reservoir in the recent years. Furthermore, the industrial sectors annually demand 152 MCM of water, while the agricultural sector is rated as the major water-demanding sector in the basin, consuming over 80% of the whole water consumed in the basin. However, the major portion of the agricultural water demands is supplied by the groundwater resources of the Basin and the remainder of the demand is supplied by the surface water stored in the reservoir, amounting to 1016 MCM. Despite the real portion of the reservoir supplying the water demand of the whole basin is not clear, this portion could be approximated according to the annual reservoir water releases in a historical period, as a method to estimate the water which should be supplied by the reservoir in different years.

In this paper, a 13-year long-term period beginning from the 2020–2021 water year, ending in the 2032–2033 water year is adopted as the near-future planning period. For estimating the water demands in each year, the 13-year historical period beginning from 2006, ending in 2019, is considered and the years included in this period are divided into two wet and normal years. Since the historical data suggest 1300 MCM is a proper figure for the water release to meet the normal water demands properly and sufficiently in the basin, the year with this volume of the water released from the reservoir is assumed as the normal year. The years in which more than 1300 MCM of water is released from the reservoir are placed in the wet division and the other years are put inside the normal division. Consequently, the years of the near-future planning period when the inflow to the reservoir is more than the inflow volumes to the reservoir in the assumed normal year were deemed to be the wet years, and the years with less than the normal inflow to the reservoir were assumed to be the normal years. The water released from the reservoir in the wet water year 2006–2007 is amounting to 1680 MCM, while the released water in the normal water year 2011–2012 is recorded to be 1110 MCM. As a result, the total water demands to be supplied by the reservoir is considered to be 1600 MCM, and 1100 MCM for the wet and normal water years, respectively. Then, the portion of each month of the wet and normal years is estimated by dividing the monthly water demands of these years into their total demand. Thereafter, the near-future 13-year planning period is also divided into two wet and normal water years, as indicated in

Table 2. The wet and normal water years of the near-future 13-year planning period and the corresponding inflows considering the inter-basin water transfer and the upstream withdrawals under the RCP4.5 scenario.

Data/Year	2020–2021	2021–2022	2022–2023	2023–2024	2024–2025	2025–2026	2026–2027	2027–2028	2028–2029	2029–2030	2030–2031	2031–2032	2032–2033
Inflow (MCM)	1327.8	1699.8	1489.9	1597.6	1746.4	1295.1	1227.7	1514.3	1403.1	1541.2	1142.5	1775.2	1298.2
wet year (W) Or normal year (N)	N	W	N	W	W	N	N	N	N	W	N	W	N

. Then, the annual domestic water is adopted to increase by 100 MCM, while the annual environmental water requirements are set to be increased by 100 MCM to cover the damages the environment suffered during the severe droughts in the past few years. Meanwhile, no increase is considered for the other water demand types assuming the water resources management plans would go well in the near-future period in the agricultural and industrial sectors. Accordingly, water demand of the wet and normal years of the near-future planning period is estimated to be increased by 200 MCM, resulting in 1800 MCM and 1300 MCM, respectively. The monthly water demands in these years are also determined considering the portion of each month from the total water demand of either of the wet or normal years of the historical period (

Table 3. Monthly water demands of the wet and normal water years assumed for the near-future planning period.

Year/Month	Oct	Nov	Dec	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep
Wet (MCM)	82.9	94.53	98.77	50	50	55.28	144.37	213.03	264.15	254.24	274.41	249.75
Normal (MCM)	77	145.4	153	55.37	52.9	78.44	149.25	179.57	166.97	77.92	83.68	80.5

).

Moreover, according to this point that in each year the reservoir must supply the domestic (drinking and sanitary) demands at an amount of 400 MCM and must also supply the demand of the industry by 40 MCM, an amount of 40 MCM is considered to be the minimum monthly water demand which must be allocated by the reservoir. This figure is increased to 50 MCM to be assigned to the near-future planning period.

It is worth mentioning that in this research, not that much drought condition was considered for the years coming in the near-future planning period. The main reason for this decision was that if the water demands, especially the agricultural water demands that are to be supplied by the reservoir was set to be at a very low level, then there would be a very large portion of the water demands to be supplied by the groundwater resources. Having no clear vision on the potential of the groundwater resources to sufficiently supply the large agricultural water demands as well as considering the risk of groundwater over-exploitation that might lead to making the vital groundwater resources unsustainable to further use in the long-term future period and to avoid considering demand management plans in the case of the groundwater overuse, it is decided not to take the drought condition into account as a climatic condition in the planning period, to impede the problem from being more complicated to solve.

The evaporation from the reservoir is rated as the main factor accounting for the water loss from the reservoir surface. In this paper, the Torrent White method is adopted to approximate the evaporation through holding relation between the evaporation and the mean temperature, as can be seen in Equations (28)–(31) [58].

$${\mathit{ET}}_p=16 \times (\frac{{10\times t}_n}{I}{)}^a$$

(28)

$$I={\displaystyle \sum }_{n=1}^{12}{\left(\frac{t_n}{5}\right)}^{1.514}={\displaystyle \sum }_{n=1}^{12}{i}_n$$

(29)

$$i_n=(\frac{t_n}{5}{)}^{1.514}$$

(30)

$$a=(0.675\times I^3-77.1\times I^2{+17920\times I+492390)\times 10}^{-6}$$

(31)

where ${\mathit{ET}}_p$ denotes the monthly potential evapotranspiration (mm); I is the annual heat index; $i_n$ is the monthly heat index; and $t_n$ represents the mean monthly temperature (°C). Equation (31) is prescribed to be used only when ${0 < t}_n \le 26.5$. When $t_n \le 0$, $E{T}_p$ is assumed to be zero, and when $t_n > 26.5$, the ${\mathit{ET}}_p$ is calculated based on the different values of $t_n$, as indicated in

Table 4. Monthly potential evapotranspiration for the months (mm) in which $t_n>26.5$ (°C).

$${\mathit{t}}_{\mathit{n}}$$	27	28	29	30	31	32	33	34	35	36	37
$${\mathit{ET}}_{\mathit{p}}$$	139	147	155	162	168	173	177	180	183	184	185

. Since these equations only calculate the evapotranspiration for the 30-day months and 12 h sunshine per day, the obtained evapotranspiration value is required to be modified for the other months and/or the other latitudes the focused study area is located in.

3.2. Reservoir Operation Optimization Model

The main goal of the optimization model solved in this paper is to maximize the sustainability index (SI) of the Zayandehrud Reservoir operation subject to several constraints mainly aimed at maintaining the reservoir sustainability at a desirable level to guarantee the water allocation from the reservoir during the planning period and even afterward. The planning period is a 13-year near-future period, beginning from the water year 2020–2021, ending in the water year 2032–2033. The sustainability index is assumed to be calculated based on Equation (33), as recommended by Sandoval-Solis et al. [37]. The performance criteria of a reservoir or any water resource system contribute to the sustainability index. These performance criteria mainly include reliability (Rel), resilience (Res), and vulnerability (Vul). However, in this paper, these criteria are calculated based on fuzzy logic concepts. Safavi and Golmohammadi [41] expressed the effectiveness of using the fuzzy-based performance criteria in holding a direct and reliable relation between the uncertainties in the water demand satisfaction and the water deficit in each time step to make a water resource system performance assessment more tangible and more reliable for the decision-makers and also the stakeholders. A fuzzy membership function can represent a degree of fuzziness of a fuzzificable variable. In this paper, the water deficit/surplus is assumed as the fuzzificable variable which is assigned the bell-shaped fuzzy membership function. This type of membership function has three tunable parameters of a, b, and c, that are set based on the experts’ points of view in this paper, as recommended in [41]. The full water deficit/surplus could be taken as the failure of a water resource system. This is why the membership degree of 0 must be possible to be assigned to this case. On the contrary, no water deficit/surplus may be taken as the full success of a water resource system and as a result, the membership degree of 1 must be assigned to this case. The bell-shaped membership function is asymptotic to 0 and also receives 1 in its range of the fuzzy membership function values, making it eligible to represent the water deficit/surplus as the fuzzificable variable to which a membership degree must be assigned. Furthermore, the fact that the bell-shaped function is asymptotic to 0 can impede the sustainability index to obtain the absolute 0 in its range of values. On the other hand, having three tunable parameters can make the bell-shaped function more flexible as compared to the other types of the membership functions such as the Gaussian membership function which only has two tunable parameters. The detailed formulation of the optimization model of the reservoir operation is as follows:

$$\mathit{Minimize} Z=\left(100-\mathit{SI}\right){\times (1+\mathit{Pen}}_1{+\mathit{Pen}}_2{+\mathit{Pen}}_3)$$

(32)

Subject to:

$$\mathit{SI}=\sqrt[3]{\mathit{Rel}\times \mathit{Res}\times \left(100-\mathit{Vul}\right) }$$

(33)

$$\begin{gathered} {\mu }_i^n=\frac{1}{1+{\left|\frac{x\times 100-c}{a}\right|}^{2b}} ; x=\mathit{Min}\left(\frac{R_i^{n }{+\mathit{SP}}_i^{n }}{{\mathit{dem}}_i^{n }},1\right) ; \\ a=30; b=2; c=100; \mathrm{for} i=1, 2,\dots , 12 , n=1, 2,\dots , 13 \end{gathered}$$

(34)

$${\mathrm{If} \mu }_i^n{<\mu }_{i+1}^{n }{ \mathrm{Then} W}_i^{n }{=\mu }_{i+1}^{n }{-\mu }_i^{n }{ \mathrm{Otherwise} W}_i^{n }=0 , \mathrm{for} i=1, 2,\dots , 11; n=1, 2,\dots , 13$$

(35)

$${\mathrm{If} \mu }_{12}^n{<\mu }_1^{n+1}{ \mathrm{Then} W}_i^{n }{=\mu }_1^{n+1 }{-\mu }_{12}^{n }{ \mathrm{Otherwise} W}_i^{n }=0; \mathrm{for} n=1, 2, \dots , 13$$

(36)

$$\mathit{Rel}=\frac{{{\displaystyle \sum }}_{n=1}^{13}{{\displaystyle \sum }}_{i=1}^{12}{ \mu }_i^{n }}{12\times 13}\times 100; 0\le \mathit{Rel}\le 100$$

(37)

$$\mathit{Res}=\frac{{{\displaystyle \sum }}_{n=1}^{13}{{\displaystyle \sum }}_{i=1}^{12}{ W}_i^{n }}{{{\displaystyle \sum }}_{n=1}^{13}{{\displaystyle \sum }}_{i=1}^{12}{(1- \mu }_i^{n })}\times 100; 0\le \mathit{Res}\le 100$$

(38)

$$\mathit{Vul}=\frac{{{\displaystyle \sum }}_{n=1}^{13}{{\displaystyle \sum }}_{i=1}^{12}{(1- \mu }_i^{n })}{12\times 13}\times 100; 0\le \mathit{Vul}\le 100$$

(39)

The spilled water in the ith month of the nth water year is denoted by ${\mathit{SP}}_i^{n }$ in Equation (34) and is calculated based on the following equations:

$$\begin{gathered} \mathrm{If} \left(S_0^1{+I}_1^1{-R}_1^1{-\mathit{loss}}_1^1\right){>S}_{\mathit{max} }{ \mathrm{Then} \mathit{SP}}_1^1{=S}_0^1{+I}_1^1{-R}_1^1{-\mathit{loss}}_1^1{-S}_{\mathit{max} } \\ {\mathrm{Otherwise} \mathit{SP}}_1^1=0 \end{gathered}$$

(40)

$$\begin{gathered} \mathrm{If} \left(S_i^n{+I}_{i+1}^n{-R}_{i+1}^n{-\mathit{loss}}_{i+1}^n\right){>S}_{\mathit{max} }{ \mathrm{Then} \mathrm{SP}}_{i+1}^n{=S}_i^n{+I}_{i+1}^n{-R}_{i+1}^n{-\mathit{loss}}_{i+1}^n{-S}_{\mathit{max} } \\ {\mathrm{Otherwise} \mathit{SP}}_{i+1}^n=0; \mathrm{for} i=1, 2,\dots , 11; n=1, 2,\dots , 13 \end{gathered}$$

(41)

$$\begin{gathered} \mathrm{If} \left(S_{12}^n{+I}_1^{n+1}{-R}_1^{n+1}{-\mathit{loss}}_1^{n+1}\right){>S}_{\mathit{max} }{ \mathrm{Then} \mathit{SP}}_1^{n+1}{=S}_{12}^n{+I}_1^{n+1}{-R}_1^{n+1}{-\mathit{loss}}_1^{n+1}{- S}_{\mathit{max} } \\ {\mathrm{Otherwise} \mathit{SP}}_1^{n+1}=0; \mathrm{for} n=1, 2,\dots ,13 \end{gathered}$$

(42)

$$\begin{gathered} { Ā }_i^n=\left(4.22\times {10}^{-17}\right){\times (S}_i^n{)}^6+\left(2.11\times {10}^{-13}\right){\times (S}_i^n{)}^5 \\ +\left(-4.25\times {10}^{-10}\right){\times (S}_i^n{)}^4+\left(4.33\times {10}^{-7}\right){\times (S}_i^n{)}^3+\left(-2.30\times {10}^{-4}\right){\times (S}_i^n{)}^2 \\ +\left(0.0916\right){\times (S}_i^n)+\left(1.087\right); \\ \mathrm{for} i=1, 2,\dots , 12; \mathrm{for} n=1, 2,\dots , 13 \end{gathered}$$

(43)

$${\mathit{loss}}_i^n{=(\mathit{evapo}}_i^n{-p}_i^n{)\times Ā }_i^n{\times 10}^{-3}; \mathrm{for} i=1, 2,\dots 12; n=1, 2,\dots ,13$$

(44)

The reservoir storage at each time step is calculated as follows:

$$S_1^1{=S}_0^1{+I}_1^1{-R}_1^1{-\mathit{loss}}_1^1{-\mathit{SP}}_1^1$$

(45)

$$S_{i+1}^n{=S}_i^n{+I}_{i+1}^n{-R}_{i+1}^n{-\mathit{loss}}_{i+1}^n{-\mathit{SP}}_{i+1}^n; \mathrm{for} i=1,\dots 11; n=1,2,\dots , 13$$

(46)

$$S_1^{n+1}{=S}_{12}^n{+I}_1^{n+1}{-R}_1^{n+1}{-\mathit{loss}}_1^{n+1}{-\mathit{SP}}_1^{n+1}; \mathrm{for} n=1, 2,\dots , 13$$

(47)

The bound constraints and the penalty function constraints of the model are formulated as follows:

$${\mathit{dem}}_{i,\mathit{min}}^{ } \le { R}_i^n\le { \mathit{dem}}_i^n; \mathrm{for} i=1, 2,\dots , 12 ; n=1, 2,\dots , 13$$

(48)

$$S_{\mathit{dead}} \le { S}_i^n\le S_{\mathit{max}}; \mathrm{for} i=1, 2,\dots , 12 ; n=1, 2,\dots , 13$$

(49)

$$S_i^n{=S}_{\mathit{opt} }; \mathrm{for} i=12; n=1, 2,\dots , 13$$

(50)

$$\begin{gathered} {\mathit{Pen}}_1{=C}_1\times {\displaystyle \sum }_{t=1}^T\left[\mathit{max}(\mathit{sgn}\left(S_{\mathit{dead}}{-S}_t\right),0)\times \left(S_{\mathit{dead}}{-S}_t\right)/S_{\mathit{dead}}\right]; \\ \mathrm{for} t=1, 2,\dots , 157 \end{gathered}$$

(51)

$$\begin{gathered} {\mathit{Pen}}_2{=C}_2\times {\displaystyle \sum }_{t=1}^T{\left[\left(S_{\mathit{opt}}{-S}_t\right)/S_{\mathit{max}}\right]}^2; \\ \mathrm{for} t=13, 25, 37, 49, 61, 73, 85, 97, 109, 121, 133, 145, 157 \end{gathered}$$

(52)

$${\mathit{Pen}}_3{=C}_3\times {\displaystyle \sum }_{t=1}^T\left[\left({\mathit{SP}}_i^n\right)/{\mathit{dem}}_i^n\right]; \mathrm{for} t=1, 2,\dots , 157$$

(53)

The parameters of the optimization model are defined as follows:

SI = Sustainability index.

${\mu }_i^n=$ Fuzzy membership degree of the reservoir performance in fuzzy demand satisfaction set in the ith month of the nth year.

$R_i^{n }=$ Water volume released from the reservoir in the ith month of the nth year (MCM).

${\mathit{SP}}_i^{n }=$ Water volume spilled from the reservoir in the ith month of the nth year (MCM).

${\mathit{dem}}_i^{n }=$ Water volume demanded downstream of the reservoir in the ith month of the nth year (MCM).

$S_i^n=$ Water volume stored in the reservoir in the ith month of the nth year (MCM).

$I_{i+1}^n=$ Inflow volume to the reservoir in the (i+1)th month of the nth year (MCM).

$R_{i+1}^n=$ Water volume releases from the reservoir in the (i+1)th month of the nth year (MCM).

${\mathit{loss}}_{i+1}^n=$ Water loss from the reservoir in the (i+1)th month of the nth year (MCM).

$S_{\mathit{max} }=$ Reservoir storage capacity set to be 1470 MCM.

${\mathit{SP}}_{i+1}^n=$ Water volume spilled from the reservoir in the (i+1)th month of the nth year (MCM).

${ Ā }_i^n=$ Area of the reservoir lake surface in the ith month of the nth year (km²).

${\mathit{evapo}}_i^n=$ Height of the stored water evaporated from the reservoir surface in the ith month of the nth year (mm).

$p_i^n=$ Height of the precipitation falling on the reservoir in the ith month of the nth year (mm).

${\mathit{dem}}_{i,\mathit{min}}^{ }=$ Minimum water demand required to be met in the ith month of the year set to 50 MCM.

$S_{\mathit{dead}}=$ Dead reservoir storage volume set to 150 MCM.

$S_{\mathit{opt}}=$ Optimal reservoir storage volume required to be kept in the reservoir at the final of each water year for preserving the sustainability of the reservoir, set to be 450 MCM.

sgn = Sign function in mathematics.

${\mathit{Pen}}_1=$ First penalty function.

${\mathit{Pen}}_2=$ Second penalty function.

${\mathit{Pen}}_3=$ Third penalty function.

$C_1=$ Penalty coefficient of the first penalty function.

$C_2=$ Penalty coefficient of the second penalty function.

$C_3=$ Penalty coefficient of the third penalty function.

i = Indicator of the month of each water year of the planning period.

n = Indicator of each water year of the planning period.

t = Indicator of the month included in the total planning period.

Z = Objective function with the penalty functions added as the constraints to be minimized.

4. Results and Discussion

In this paper, the time series of all the input data necessary to launch the optimization process of the reservoir operation including the precipitation, temperature, and the inflow to the reservoir are extracted from [48,49,50], and then the GSA and the proposed AAGSA optimization models are implemented in the future period beginning from 2020–2021, ending in 2032–2033, under the climate change RCP4.5 and RCP8.5 scenarios. The algorithms in both scenarios are run for 30 independent runs. The minimum fitness function value obtained in the best run of either of the examined algorithms is adopted as the final result of the optimization process, resulting in the best decision variables (water releases) and state variables (water storages) for the case study.

As the results suggest, for RPC4.5, the best fitness function value given by the proposed AAGSA is 56 times (98.2%) less than that offered by the GSA, while for RPC8.5, this value is 44 times (97.7%) less than that evaluated by the GSA. Furthermore, the runtime of the AAGSA is nearly 5 times (80%) less than that of the GSA in both scenarios, suggesting the best overall performance of the AAGSA compared to the GSA in terms of both the quality of the results and the computational cost. Figure 3 and Figure 4 show the convergence curves plotted upon implementing the AAGSA against those plotted when running the GSA for RCP4.5 and RCP8.5, respectively. The detailed parameter settings of the algorithms and the final results obtained upon running both the GSA and AAGSA algorithms for RCP4.5 and RCP8.5 are presented in

Table 5. Parameter setting and the final results of the GSA and AAGSA when applied to the reservoir operation problem under RCP4.5 and RCP8.5.

Parameters/Results	Characteristics	RCP4.5		RCP8.5
		AAGSA	GSA	AAGSA	GSA
Parameters	Total number of iterations	1000	1000	1000	1000
	Number of search agents	50	50	50	50
	Total number of runs	30	30	30	30
	${\mathit{final}}_{\mathit{per}}$	2	2	2	2
	$G_0$	100	100	100	100
	α	20	20	20	20
	$C_{\mathit{random}}$	2	-	2	-
	${\mathit{mut}}_{\mathit{max}}$	0.9	-	0.9	-
	${\mathit{mut}}_{\mathit{min}}$	0.1	-	0.1	-
	${\mathit{power}}_{\mathit{max}}$	4	-	4	-
	${\mathit{power}}_{\mathit{min}}$	1	-	1	-
Results	Run time (seconds)	1463.37	7371.63	1376.49	6216.59
	First penalty function value	0	0	0	0
	Second penalty function value	0.8736	2.5027	3.0476	4.052
	Third penalty function value	3.0917	16.7218	32.1338	55.3013
	Minimum fitness function value	7.3083	408.551	19.6172	864.8347
	Average fitness function value	18.4301	463.2886	68.7664	1040.2989
	Standard deviation of fitness function v Values over the total runs	7.9176	40.5558	23.2623	99.0883
	Median fitness function values	16.7354	456.5743	62.4131	1045.0276
	Sustainability index (%)	98.53	79.79	99.46	85.67

.

The detailed results of applying both algorithms to the reservoir operation problem for RCP4.5 and RCP8.5 are presented in

Table 6. Detailed results showing the performance criteria of the GSA and AAGSA in solving the reservoir operation problem under RCP4.5.

Criteria	Algorithm	2020–2021	2021–2022	2022–2023	2023–2024	2024–2025	2025–2026	2026–2027	2027–2028	2028–2029	2029–2030	2030–2031	2031–2032	2032–2033
Demand met (%)	AAGSA	91.65	92.31	92.32	94.2	94.71	93.7	96	95.2	96.29	97.21	94.13	92.77	92.9
	GSA	81.49	83.56	87.05	87.91	84.56	83.76	85.77	83.34	88.7	85.38	85.2	79.35	81.57
Vulnerability	AAGSA	2.71	2.49	2.3	2.25	3.09	1.52	0.32	0.29	0.18	0.35	0.94	5.81	2.04
	GSA	20.96	20.06	12.78	14.68	18.61	15.62	16.6	20.23	9.94	20.87	14.71	32.95	20.53
Resilience	AAGSA	99.26	98.78	98.91	99.1	99.92	99.12	98.33	95.45	97.72	99.87	99.62	99.94	98.89
	GSA	63.72	75.48	94.54	74.83	86.38	83.1	94.38	63.82	86.5	75.32	92.03	62.13	69.61
Reliability	AAGSA	97.29	97.51	97.7	97.75	96.91	98.48	99.68	99.71	99.82	99.65	99.06	94.19	97.96
	GSA	79.04	79.94	87.22	85.32	81.39	84.38	83.4	79.77	90.06	79.13	85.29	67.05	79.47
SI	AAGSA	97.94	97.93	98.1	98.2	97.9	98.69	99.23	98.27	99.11	99.72	99.25	96.07	98.27
	GSA	73.57	78.43	89.59	81.67	83.02	83.95	86.91	74.05	88.86	77.84	87.48	65.37	76.03

and

Table 7. Detailed results showing the performance criteria of the GSA and AAGSA in solving the reservoir operation problem under RCP8.5.

Criteria	Algorithm	2020–2021	2021–2022	2022–2023	2023–2024	2024–2025	2025–2026	2026–2027	2027–2028	2028–2029	2029–2030	2030–2031	2031–2032	2032–2033
Demand met (%)	AAGSA	92.78	95.42	99.31	99.81	96.83	95.15	96.65	95.96	99.70	96.78	96.02	97.42	95.97
	GSA	84.18	85.93	88.86	90.24	89.10	92.22	94.74	89.96	86.73	86.83	86.50	87.45	87.64
Vulnerability	AAGSA	1.94	1.12	0.00	0.00	1.48	0.82	0.61	0.64	0.00	0.93	0.58	0.56	0.37
	GSA	17.40	18.79	11.86	13.99	14.62	5.90	2.72	10.51	14.88	17.78	13.50	15.66	11.50
Resilience	AAGSA	99.72	99.71	99.62	100.00	99.89	99.66	99.79	99.76	100.00	99.82	99.80	99.78	99.72
	GSA	94.26	92.46	95.59	70.70	50.64	95.84	64.02	90.30	97.95	74.39	83.64	78.44	85.57
Reliability	AAGSA	98.06	98.88	100.00	100.00	98.52	99.18	99.39	99.36	100.00	99.07	99.42	99.44	99.63
	GSA	82.60	81.21	88.14	86.01	85.38	94.10	97.28	89.49	85.12	82.22	86.50	84.34	88.50
SI	AAGSA	98.61	99.16	99.87	100.00	98.98	99.34	99.52	99.49	100.00	99.32	99.55	99.55	99.66
	GSA	86.32	84.80	90.56	80.57	71.74	94.68	84.62	89.76	89.20	79.52	85.54	82.32	87.51

, respectively. As

Table 8. Detailed results showing the overall performance criteria of the GSA and AAGSA algorithms in solving the reservoir operation problem over the whole planning period under RCP4.5 and RCP8.5.

Climate Change Scenario	Algorithm	Sustainability Index	Reliability	Resilience	Vulnerability	Water Demand Met (%)
RCP4.5	AAGSA	98.53	98.13	99.32	1.87	94.10
	GSA	79.80	81.65	76.22	18.35	84.43
RCP8.5	AAGSA	99.46	99.30	99.77	0.70	96.75
	GSA	85.67	86.99	83.09	13.01	88.49

suggests, the average performance criteria of the AAGSA over the whole planning period is significantly improved as compared to those of the original GSA for both climate change scenarios. In detail, when employing the proposed AAGSA against the original GSA for solving the problem under the RCP4.5 scenario, the vulnerability index is decreased by 89.8%, resilience index is increased by 30%, reliability index is increased by 20%, and the final sustainability index is increased by 23.5%. While implementing the AAGSA against the GSA to solve the problem under the RCP8.5 scenario reveals that the vulnerability index is decreased by 94.6%, resilience index is increased by 20%, reliability index is increased by 14%, and the final sustainability index is increased by 16%. Moreover, the average water demand percentages met calculated by the proposed AAGSA are increased by 11.5% and 9.3% as compared to those obtained by the original GSA when solving the reservoir operation problem under the RCP4.5 and RCP8.5 scenarios, respectively.

As can be seen in Table 8, the performance criteria of the AAGSA and GSA under RCP8.5 are slightly improved compared to RCP4.5. In addition, the improvement obtained is more significant in the performance of the GSA than that of the AAGSA. The main reason for this improvement is hidden in the natural characteristics of the RCP8.5 against the RCP4.5. Under the RCP8.5, the averaged inflow to the reservoir in the future planning period is predicted to be increased by 24.96% and the average temperature is estimated to increase by 2.3% [48,49,50]. The highly increased inflow to the reservoir might have led to a higher percentage in meeting the water demands in both algorithms alongside the higher amounts for reservoir storage. While the conditions are eased for either of the algorithms to solve such a problem, the major weakness of the GSA in the exploration of the search space is highly downgraded as the optimum is not that extreme in the search space and gets closer to the center of the problem’s domain, entailing any algorithm to have more effective exploitation capability than the higher exploration skill. Whereas the performance of the AAGSA remains strong and efficient when the scenario is switched from the RCP4.5 to the RCP8.5, as the proposed AAGSA has high competence in both the exploration and exploitation that allows this algorithm to easily solve the newly emerging problem, as well. As the results presented in Table 8 suggest, the SI for AAGSA and GSA under RCP8.5 is increased by 0.94% and 7.4%, respectively, as compared to those obtained for these algorithms under RCP4.5. Moreover, the water demand percentage met obtained under RCP8.5 is increased by 2.4% and 4.8% as compared to those yielded under RCP4.5, as calculated by the AAGSA and GSA, respectively.

Figure 5, Figure 6, Figure 7 and Figure 8 show the monthly volume of the water released, the monthly volume of the downstream water demand, and the monthly volume of the inflow to the reservoir over the whole planning period resulting from the GSA and AAGSA algorithms implemented under RCP4.5 and RCP8.5. As these figures suggest, the AAGSA can better supply the demands over the planning period. The main reasons for the AAGSA to more appropriately meet the water demands can be divided into three main points: (1) the reservoir storage is at a desirably high level to meet the demands, especially in the dry periods; (2) the inflow to the reservoir is properly high within the periods with high rates of the precipitation and the resulting runoff volumes; and (3) the water demands are essentially at low levels at the months whose water demands are desirably met. Among these three reasons, the first one is the factor related to the performance of the optimization models as the experiments accomplished in this paper suggest the superiority of the AAGSA to handle the optimization problems of any type and any condition. The second and the third reasons are related to the conditions experienced by the problem in the focused case study and are both consequently assumed constant for the GSA and the proposed AAGSA. As a result, an algorithm could have better performance than the other one when getting more adapted to these conditions when operating the reservoir. As the detailed results obtained under RCP4.5 indicate, the AAGSA and GSA could meet more than 95% of the whole monthly water demands in 100 out of 156 (64%), and 43 out of 156 (27.6%) of the planning period, respectively. Under RCP8.5, these figures are turned into 78.85% and 42.95%, for the AAGSA and GSA, respectively. Meanwhile, both algorithms can meet the water demands in nearly 95% of the months having been assigned the least demands under both climate change scenarios. Moreover, the AAGSA could meet the demands by more than 85% in all the months after the months May and June of all years of the planning period when the reservoir storage is at the maximum level in both scenarios, while the GSA could meet the demands of such months by more than 85% at only 11 out of 26 (42.3%) of these months, suggesting the much worse performance of the GSA to make the water demands sustainable to be met over a critical future planning period under a climate change scenario; however, this matter occurs in 21 out of 26 (80.8%) of these months under RCP8.5, which is due to the GSA’s better performance under this scenario as previously mentioned.

As Figure 9 and Figure 10 display, when the reservoir operation problem is solved by the GSA, the reservoir tends to hold its storage volume at much higher levels than that offered by the AAGSA. This performance of the GSA can make the releases from the reservoir less than those presented by the AAGSA, and the spilled water volume much more than that offered by the AAGSA. This behavior of the GSA to handle the reservoir operation problem impedes this algorithm to generate the favorite results, at least as compared to the proposed AAGSA, as the goals of the optimization model are (1) maximizing the sustainability of the water demands to be met, and (2) holding the reservoir storage at a well-balanced level, such that the storage level is neither too low to make the water demands in the dry months of the planning period not be desirably met, nor too high to frequently lead the water to be spilled from the reservoir, as the monthly spilled water volume is restricted in the formulation of the optimization model. The latter point is realized to not be properly met by the GSA, as this algorithm makes the water storage of the reservoir be spilled at 22 out of 156 months (14%) of the planning period under RCP4.5 and 53 out of 156 months (34%) under RCP8.5, while the proposed AAGSA only makes the water be spilled at 4 out of the 156 months (2.5%) of the total planning period under RCP4.5 and 35 out of 156 months (22.4%) under RCP8.5. As previously mentioned, the reason for higher levels of reservoir storage volumes under RCP8.5 compared to those under RCP4.5 is referred to highly increased inflow volumes in the former scenario [48,49,50].

5. Conclusions

The main goal of this study was to maximize the sustainability index of a reservoir to meet the downstream water demands. The optimization model was also subject to several constraints the major of which is to keep the reservoir storage sustainable during the whole planning period to provide better conditions for the reservoir to satisfy the water demands as much as possible. A future period beginning from 2020–2021, ending in 2032–2033, was assumed as a 13-year long-term planning period whose data necessary to simulate its climatic conditions were obtained as the outputs of the RCP4.5 and RCP8.5 climate change scenarios. The RCP4.5 was chosen as the climate change scenario mainly because of two main points: (1) the RCP climate change scenarios are not so different in the outputs and the data they project in the near-future periods such as the period adopted to be the planning period for optimizing the reservoir operation in this paper, and (2) the RCP4.5 was adopted as a medium-standing scenario and thus, a representative scenario for all the other RCP scenarios. As a result, the examination of the effects of the other scenarios was avoided in the near-future planning period. After the results of solving the reservoir operation problem were revealed under the RCP4.5 scenario, the two optimization models were also applied to the problem under the RCP8.5 scenario, to benchmark the possible differences between these two scenarios and their effect on the way to operate the reservoir during the near-future planning period. The examination of the RCP8.5 could also help how one can rely on the RCP4.5 as a representative for all the other climate change scenarios.

The classic sustainability index involves Boolean logic in its calculations. The Boolean logic rates the performance of a focused water resource system as a success only when fully supplying the water demands, and only considers it as a failure when not managing to fully supply the water demands by that focused water resource system. This behavior of the Boolean logic cannot provide an accurate, tangible, and plausible view for the decision-makers and stakeholders when the reservoir operates. To avoid these problems, the fuzzy logic is used to calculate the components of the sustainability index in this paper, as the fuzzy logic partly validates the performance of a water resource system when supplying the water demands at any level possible for that system. The Zayandehrud Reservoir was adopted as the case study of this paper. This reservoir has experienced severe climatic changes, droughts, and mismanagements during recent years and thus it is necessary to adapt the reservoir to the current trend of climatic conditions which can be projected by the climate change scenarios, helping the operating policies of such a reservoir be optimally updated. An adaptive accelerated gravitational search algorithm (AAGSA) is developed as a new variant of the robust gravitational search algorithm (GSA) to solve the optimization model built up in this paper. The AAGSA has two main advantages against the GSA, discriminating it from the GSA and making it a highly effective meta-heuristic algorithm, especially to solve the large-scale optimization problems: (1) the exploration phase of the AAGSA is significantly accelerated to avoid the algorithm to be trapped in local optima and to entirely search the decision space to detect all the promising regions in this space before the exploitation process is started, and (2) an effective balance is held in the AAGSA between the exploration and exploitation phases of the optimization process through de-emphasizing the masses (fitness) of the search agents at the early iterations while emphasizing the masses and de-emphasizing the distances between the agents up to the final iterations. The AAGSA demonstrated higher-quality results while reducing the runtime by 80%, as compared to the original GSA. The GSA presented higher storage volumes for the reservoir causing the reservoir to spill the water, whereas this behavior is rated undesirable during a reservoir operation period and thus, the reservoir was restricted to unduly spill the water as formulated in the optimization model. Furthermore, under RCP4.5 and RCP8.5, the proposed AAGSA made the reservoir averagely meet the water demands by 94% and 96.75% (11.5% and 9.3% better than that offered by the GSA) and present the total sustainability index to be 98.5% and 99.46% (23.5% and 16% better than that suggested by the GSA). Finally, the AAGSA improved the fitness function of the optimization by 98.2% compared to the GSA under RCP4.5, and by 97.7% under RCP8.5 compared to the GSA, demonstrating the significant superiority of the proposed AAGSA to the original GSA in solving such a high-dimensional and highly constrained practical engineering problem.