This work introduces the symbiotic organisms search (SOS) evolutionary algorithm to the optimization of reservoir operation. Unlike the genetic algorithm (GA) and the water cycle algorithm (WCA) the SOS does not require specification of algorithmic parameters. The solution effectiveness of the GA, SOS, and WCA was assessed with a single-reservoir and a multi-reservoir optimization problem. The SOS proved superior to the GA and the WCA in optimizing the objective functions of the two reservoir systems. In the single reservoir problem, with global optimum value of 1.213, the SOS, GA, and WCA determined 1.240, 1.535, and 1.262 as the optimal solutions, respectively. The superiority of SOS was also verified in a hypothetical four-reservoir optimization problem. In this case, the GA, WCA, and SOS in their best performance among 10 solution runs converged to 97.46%, 99.56%, and 99.86% of the global optimal solution. Besides its better performance in approximating optima, the SOS avoided premature convergence and produced lower standard deviation about optima.

## INTRODUCTION

The unwise operation of reservoirs is the main driving-force of various water resources crises such as degrading native aquatic ecosystems (Steinschneider *et al.* 2014), water pollution (Yuan *et al.* 2015), shrinkage of lakes (Azarnivand & Banihabib 2016), and other calamities in many regions of the world. The countries within the arid regions of the world are grappling with anthropogenic and climatic driving forces which pose a burden to water, energy, and food security. For this reason, it is vital to improve water resources planning and management, which includes reservoir operation as a key component.

The techniques used for obtaining optimal operation of reservoir systems can be divided into classical methods and evolutionary algorithms (EAs). Schardong & Simonovic (2015) warned about the curse of dimensionality that plagues classical methods such as dynamic programming (DP) and stochastic dynamic programming (SDP). Linear programming (LP) requires objective functions and constraints that are linear on the decision variables, this being a common limitation for practical modeling of real reservoir systems. Nonlinear programming (NLP) is heavily influenced by the choice of initial conditions, and most searches may lead to local optimal solutions instead of the global optima. On the other hand, EAs have amply demonstrated their ability of finding near-global solutions of complex optimization problems (Farhangi *et al.* 2012; Zufferey 2012; Ashofteh *et al.* 2013a, 2013b, 2015a, 2015b, 2015c; Bozorg-Haddad *et al.* 2013, 2014, 2015a, 2015b; Fallah-Mehdipour 2013a, 2013b, 2013c, 2014; Orouji *et al.* 2013, 2014a, 2014b; Shokri *et al.* 2013; 2014; Soltanjalili *et al.* 2013; Ahmadi *et al.* 2014, 2015; Beygi *et al.* 2014; Bolouri-Yazdeli *et al.* 2014; Jahandideh-Tehrani *et al.* 2015). There is a large volume of published papers dealing with EAs applied to reservoir operation problems (Ching *et al.* 2003; Kumphon 2013; Afshar *et al.* 2015; Bashiri-Atrabi *et al.* 2015; Li *et al.* 2015; Rampazzo *et al.* 2015; Amirkhani *et al.* 2016; Bozorg-Haddad *et al.* 2016a, 2016b).

An EA such as the standard genetic algorithm (GA) can be applied to solving many types of optimization problems. The GA mimics natural selection mechanisms and works on the basis of populations of solutions that are improved iteratively (i.e., population-by-population (Fogel 2000)). In the first computational step of the GA, an initial population is created composed of randomly generated solutions. The next population is produced to improve the objective function through an iterative process. The chromosomes (or solutions) of the current population at each step are selected to generate the next generation. The selection probability of chromosomes with superior fitness is larger than those of less fit chromosomes. Selected chromosomes generate the next population with crossover and mutation operators. Crossover generates two new chromosomes by exchanging genes between them. The mutation alters the chromosomes' genes to create diversity in their population. The process of improving chromosomal populations continues iteratively until fulfilling a termination criteria.

Based on the ‘no free-lunch’ theorem, it is impossible for an EA to optimally solve all the optimizing problems (Wolpert & Macready 1997). Thus, there are new EAs being introduced continually. Yazdi *et al.* (2016) compared three EAs for the optimal design of buildings. Hosseini-Moghari *et al.* (2015) compared a recently developed EA called imperialist competitive algorithm (ICA) against the cuckoo optimization algorithm (COA) in two separated optimization problems, in which ICA outperformed COA in solving single reservoir and multi-reservoir operation problems. The ICA also outperformed ant colony optimization (ACO) in deriving optimal operational policies of the Dez reservoir in Iran (Afshar *et al.* 2015). Li *et al.* (2015) compared seven typical heuristic algorithms with an application to a real-life multi-reservoir system in China. In that study, particle swarm optimization proved superior to other algorithms. Bashiri-Atrabi *et al.* (2015) developed an optimization model based on the comparison of the harmony search algorithm vs. the honey-bee mating optimization for reservoir operation optimization with respect to flood control in northern Iran. Akbari-Alashti *et al.* (2014) applied NLP, GA, and fixed-length gene genetic programming (FLGGP) to derive multi-reservoir real-time operation rules for a multi-reservoir system. The results indicated the superiority of FLGGP in reaching the global optimum value of NLP.

The selection of an EA to solve a specific problem remains non-trivial (Maier *et al.* 2014). EAs are targets of criticism because of the need for specifying algorithmic parameters. Parameter tuning (either implicitly or explicitly) of EAs can be complex and time-consuming (Lobo *et al.* 2007; Eiben & Smit 2011; Joan-Arinyo *et al.* 2011; Yeguas *et al.* 2014; Veček *et al.* 2016). Parameter tuning using full factorial design is computationally burdensome. Hence, practitioners often apply already-available tuning approaches or develop parameter specification methods (Joan-Arinyo *et al.* 2011; Lee *et al.* 2013; Montero *et al.* 2014; Veček *et al.* 2016). It appears practical to test an EA which does not involve parameter tuning in water resources management. Cheng & Prayogo (2014) introduced the symbiotic organisms search (SOS) algorithm to overcome the parameter-specification disadvantage. The SOS algorithm requires only the specification of the ‘maximum number of evaluations’ and the ‘population size’. The SOS is a nature-inspired optimization algorithm which simulates three different symbiosis interactions between organisms that dwell in an ecosystem.

This work focuses on the optimal hydropower production optimization of the Karun4 reservoir system in southwestern Iran. The novelty of this work consists in the introduction of the SOS algorithm, a recently developed EA which does not involve parameter tuning to water resources systems analysis. This work assesses the performance of the SOS algorithm in optimizing a single reservoir operation against the GA, the water cycle algorithm (WCA), and NLP. Moreover, the optimal operation of a hypothetical benchmark four-reservoir system is also tackled with SOS, GA, WCA, and NLP to verify the used EAs. The remainder of this paper first describes the SOS algorithm and its theoretical underpinnings. The SOS, GA, WCA, and NLP are then verified with a hypothetical four-reservoir system. Lastly, the EAs and NLP are applied to solving for the optimal optimization of the Karun4 reservoir of Iran. The performances of the SOS, GA, WCA, and NLP are compared using the solutions to the single-reservoir problems.

## MATERIAL AND METHODS

### The principles of the SOS algorithm

*X*is identified. Two randomly selected organisms along with

_{best}*X*participate in a dialectic relationship that is beneficial for both of them and is called mutualism. The following equations generate new candidate solutions: where the mutual vector (

_{best}*MV*) is the average value of

*x*and

_{i}*x*which allows the organisms to be updated simultaneously rather than separately; rand (0,1) is a vector of random numbers. In a mutualistic symbiosis between two organisms in nature, one organism might receive a great benefit while the other receives no significant benefit. This is reflected by

_{j}*BF*

_{1}and

*BF*

_{2}, which are determined randomly as either 1 or 2 ([

*BF*=

_{i}*Round*(rand (0, 1) +1];

*i*= 1 and 2) to reflect the level of benefits received from the engagement.

*x*and

_{i}*x*are then compared to the old ones. Fitter organisms are selected as solutions for the next iteration. The comparisons and selections starts with the counter 1 and ends with the counter equal to the

_{j}*population size*(

*npop*). For each

*i,*the solution

*j*is selected randomly within the current population. Then, fitter organisms participate in the next phase of interaction which is called commensalism. According to commensalism one organism receives benefits while the other remains unaffected. Similar to the previous section

*, x*is selected randomly from the ecosystem to interact with

_{j}*x*;

_{i}*x*tries to garner benefits from the relationship, yet

_{i}*x*remains neutral or unaffected. In this phase,

_{j}*x*is updated as follows, if the new fitness value outperforms the previous fitness values:

_{i}In the third phase, which entails the mutation operator of the SOS and is called parasitism, *x _{i}* and

*x*are the artificial parasite and host, respectively. In this type of symbiosis relationship, one organism benefits while another one is harmed. The trademark of the parasite vector (

_{j}*PV*) is that it competes against other randomly selected dimensions rather than its parent/creator with a range between given lower and upper bounds. In this phase, an initial

*PV*is generated by duplicating organism

*x*. Some of the decision variables from the

_{j}*PV*are modified randomly to distinguish the

*PV*from

*x*. A random number must be generated in the range of [1, the number of decision variables] to represent the total number of modified variables. A uniform random number is generated for each dimension to obtain the location of the modified variables. Lastly, a uniform distribution within the search space is required to modify the variables and provide a

_{j}*PV*for the parasitism phase. In synthesis the

*PV*attempts to replace

*x*which is selected randomly from the ecosystem. If the

_{j}*PV*outperforms

*x*it becomes part of the ecosystem, whereas if the

_{j}*PV*does not outperform

*x*it vanishes from the ecosystem. The

_{j}*PV*is created by modifying

*x*in random dimensions with random numbers instead of making small changes in

_{j}*x*. If the current

_{j}*PV*and

*x*are not the last member of the ecosystem the algorithm returns to the step that selected

_{j}*X*until reaching a specified termination criterion.

_{best}### Simulation model for reservoir operation

*t*= the index for simulation periods;

*i*= reservoir index;

*S*

_{(i,t)}and

*S*

_{(i,t+1)}= the storages of the

*ith*reservoir, respectively, at the beginning and end of period

*t*(MCM);

*Q*

_{(i,t)}= inflow volume into the

*ith*reservoir during period

*t*(MCM/s);

*M*

_{(i,t)}= the matrix of input–output connectivity between reservoirs;

*R*

_{(j,t)}= release volume from the

*jth*reservoir during period

*t*(MCM);

*Sp*

_{(j,t)}= overflow volume from the

*jth*reservoir during period

*t*(MCM);

*Loss*

_{(i,t)}= evaporation loss from the

*ith*reservoir surface during period

*t*(MCM);

*n*= number of reservoirs; and

*T*= total number of operation periods (month).

*Loss*

_{(i,t)}]: where

*Ev*

_{(i,t)}= net evaporation (evaporation minus precipitation) from the

*ith*reservoir surface during the period

*t*(Km); and reservoir areas, respectively, at the beginning and end of the period

*t*(Km

^{2}). is evaluated with the area–storage formula as follows:

*S*

_{max(i,t)}= the maximum storage in the

*ith*reservoir during period

*t*.

*R*

_{min(i,t)}and

*R*

_{max(i,t)}

*=*minimum and maximum allowable release from the

*ith*reservoir during period

*t*, respectively;

*S*

_{min(i,t)}= minimum storage value of the

*ith*reservoir at the beginning of period

*t*;

*S*

_{(i,1)}= the storage of the

*ith*reservoir in period 1;

*S*

_{(i,T+1)}= the storage of the

*ith*reservoir at the end of period

*T*.

### Case study: optimal operation of the Karun4 reservoir system

*PPC*) equals 1,000 megawatts; the minimum and maximum reservoir storages are 1,141 × 10

^{6}and 2,190 × 10

^{6}m

^{3}, respectively. The simulation–optimization model for optimal operation of the Karun4 reservoir was structured for a monthly time step during the period 1996–1997 to 2000–2001. The average volume of inflow plus evaporation depth from the surface of Karun4 reservoir during the study period is presented in Figure 3. Power production at the Karun4 reservoir is given by the following formula: where

*P*

_{(t)}= hydropower generation in period

*t*(MW);

*g*= acceleration of gravity (m/s

^{2});

*η*= efficiency of the power plant;

*Rp*

_{(t)}= release of water through the power plant in period

*t*(MCM);

*PF =*plant functional coefficient;

*Mul*

_{(t)}= 10

^{6}times of the number of seconds in period

*t*;

*H*and

_{t}*H*

_{t+1}= reservoir water level at the beginning and end of period

*t*(m), respectively;

*Tw*

_{(t)}= reservoir tail-water level which is assumed constant for all periods during period

*t*(m); and

*PPC*= power plant capacity.

*t*(see Equation (11)). The objective function is given by:

*P*

_{1}and

*P*

_{2}on storage deviations outside the feasible region are: where

*P*1 and

*P*2

_{(t)}are penalty functions on storage not meeting the carryover constraints, and being less than the minimum storage, respectively;

*K*

_{1}, and

*K*

_{2}= constants of the penalty functions.

### Verification of the algorithm with the benchmark problem: optimal operation of a four-reservoir system

*M*is a fourth-order matrix with −1 s along the diagonal denoting releases from reservoirs, and off-diagonal +1 s denoting transfers of water from one reservoir to another.

The EAs' computations for the benchmark functions, the operation of the Karun4 reservoir system, and the operation of the four-reservoir system were programmed with MATLAB (*MATLAB* 7.11.0 software). The NLP solution was obtained with Lingo based on generalized reduced gradient algorithm (Lingo 11.0 software). Multiple solutions (runs) are required to assess statistically the performance of the EAs due to their random nature. Hence, the results were obtained based on 10 runs for each applied EA.

## RESULTS AND DISCUSSION

The SOS algorithm's results for the Karun4 reservoir operation were compared with those obtained with the GA, WCA, and NLP. The NLP method was applied to evaluate the global optimal solution whereas the EAs determined the near-optimal solution. The crossover rate, mutation rate, and number of populations for GA were determined by trial and error identical to the approach by Bozorg-Haddad *et al.* (2011) and set equal to 0.60 (via two-point crossover function), 0.05 (via uniform function), and 70,000 (including population size [ecosystem size in the SOS] = 70 and number of generations =1,000), respectively. Moreover, the selection process applied the roulette wheel. Due to critique associated with the application of the roulette wheel by Bozorg-Haddad *et al.* (2011), the authors also tested other operators. The crossover rate and mutation rate were set equal to 0.70 and 0.06, respectively. The parameter values of the WCA were identical to those used by Bozorg-Haddad *et al.* (2015c). The number of objective function evaluations (NFE) for the four-reservoir problem equaled 500,000 with the EAs (the SOS algorithm, the WCA, and the GA).

### Verification of results with the benchmark problems

#### The four-reservoir system operation

*et al.*(2011) calculated the global optimal NLP solution to be 308.29. In this case, the constraints and the objective function are linear-based, which means that the NLP is truly a LP problem. Table 1 lists the results for the 500,000 objective function evaluations, from which it is deduced that the average objective-function values calculated with the SOS algorithm, WCA, and the GA equaled 306.50, 304.92, and 299.70, respectively. Thus, optimization with the SOS algorithm outperformed the WCA and the GA. (1) The superiority of SOS algorithm to the optimization methods was also true for the best performance of the three EAs. The best objective function values calculated with the SOS algorithm, WCA, and the GA were 307.85, 306.92, and 300.47, respectively. (2) The best performances of the WCA and the GA converged to 99.56% and 97.46% of the global optimum solution, respectively, while the SOS algorithm in its highest performance reached 99.86% of the global optimal solution. The convergence history of the SOS algorithm and the GA based on their average performance for the four-reservoir system is illustrated in Figure 5. Figure 5 clearly indicates the superiority of the SOS over the GA. Moreover, Figure 6(a) reveals monthly releases and storages in each reservoir, respectively, on the basis of the best optimal values from the SOS algorithm and the GA. Figure 6(a) depicts appropriate compatibility between the SOS and the NLP results. The incompatibility of optimal releases between the GA and NLP for all four reservoirs is clear from Figure 6(a). Figure 6(b) shows that the incompatibility of the optimal storage calculated with the GA and NLP for the third and fourth reservoirs is higher than that between the first and second reservoirs. A slight incompatibility is also observed between the SOS algorithm and NLP's results for the third and fourth reservoirs.

Number of run . | GA^{a}
. | WCA^{a}
. | SOS . | NLP^{b}
. |
---|---|---|---|---|

1 | 300.42 | 306.83 | 305.99 | 308.29 |

2 | 298.89 | 302.40 | 306.45 | |

3 | 300.09 | 303.65 | 307.30 | |

4 | 300.47 | 303.60 | 306.11 | |

5 | 298.46 | 302.38 | 307.85 | |

6 | 300.00 | 306.01 | 305.67 | |

7 | 299.22 | 304.05 | 306.86 | |

8 | 299.87 | 306.75 | 307.28 | |

9 | 299.20 | 306.63 | 305.47 | |

10 | 300.35 | 306.92 | 305.97 | |

Best | 300.47 | 306.92 | 307.85 | |

Worst | 298.46 | 302.38 | 305.47 | |

Average | 299.70 | 304.92 | 306.50 | |

Standard deviation | 0.7060 | 1.8863 | 0.7914 | |

Coefficient of variation | 0.0024 | 0.0062 | 0.0026 |

Number of run . | GA^{a}
. | WCA^{a}
. | SOS . | NLP^{b}
. |
---|---|---|---|---|

1 | 300.42 | 306.83 | 305.99 | 308.29 |

2 | 298.89 | 302.40 | 306.45 | |

3 | 300.09 | 303.65 | 307.30 | |

4 | 300.47 | 303.60 | 306.11 | |

5 | 298.46 | 302.38 | 307.85 | |

6 | 300.00 | 306.01 | 305.67 | |

7 | 299.22 | 304.05 | 306.86 | |

8 | 299.87 | 306.75 | 307.28 | |

9 | 299.20 | 306.63 | 305.47 | |

10 | 300.35 | 306.92 | 305.97 | |

Best | 300.47 | 306.92 | 307.85 | |

Worst | 298.46 | 302.38 | 305.47 | |

Average | 299.70 | 304.92 | 306.50 | |

Standard deviation | 0.7060 | 1.8863 | 0.7914 | |

Coefficient of variation | 0.0024 | 0.0062 | 0.0026 |

#### The Karun4 reservoir system operation

The NLP solution was equal to 1.213 (Bozorg-Haddad *et al.* 2015c). The superiority of an EA over other ones can be obtained based on the similarity of its near-optimal solution to the global optimal solution. Each run of the SOS algorithm lasted approximately 30 seconds. Table 2 demonstrates the performance of the SOS algorithm vs. GA and WCA based on 10 runs. The main results to emerge from Table 2 are as follows. (1) The GA in its best performance converged to 1.535, while the SOS algorithm reached the value 1.240 in its best performance. The SOS also outdid the WCA, which converged to 1.260 in its best run. (2) It is noteworthy that even the worst performance of the SOS algorithm was better than the best performance of the GA. (3) The coefficient of variation of the GA's solutions was almost five times greater than that of the SOS algorithm. Yet, on the basis of standard deviation and coefficient of variation, the WCAs' runs exhibited smaller variability than the SOS algorithm. The coefficients of variation of the SOS algorithm and WCA were close to zero.

Number of run . | GA^{a}
. | WCA^{a}
. | SOS . | NLP^{a}
. |
---|---|---|---|---|

1 | 1.673 | 1.289 | 1.257 | 1.213 |

2 | 1.549 | 1.269 | 1.253 | |

3 | 1.865 | 1.287 | 1.240 | |

4 | 1.752 | 1.260 | 1.291 | |

5 | 1.987 | 1.289 | 1.242 | |

6 | 1.753 | 1.285 | 1.245 | |

7 | 1.931 | 1.281 | 1.248 | |

8 | 1.57 | 1.279 | 1.314 | |

9 | 1.842 | 1.286 | 1.272 | |

10 | 1.535 | 1.262 | 1.248 | |

Best | 1.535 | 1.260 | 1.240 | |

Worst | 1.987 | 1.289 | 1.314 | |

Average | 1.746 | 1.279 | 1.261 | |

Standard deviation | 0.162 | 0.010 | 0.024 | |

Coefficient of variation | 0.093 | 0.008 | 0.019 |

Number of run . | GA^{a}
. | WCA^{a}
. | SOS . | NLP^{a}
. |
---|---|---|---|---|

1 | 1.673 | 1.289 | 1.257 | 1.213 |

2 | 1.549 | 1.269 | 1.253 | |

3 | 1.865 | 1.287 | 1.240 | |

4 | 1.752 | 1.260 | 1.291 | |

5 | 1.987 | 1.289 | 1.242 | |

6 | 1.753 | 1.285 | 1.245 | |

7 | 1.931 | 1.281 | 1.248 | |

8 | 1.57 | 1.279 | 1.314 | |

9 | 1.842 | 1.286 | 1.272 | |

10 | 1.535 | 1.262 | 1.248 | |

Best | 1.535 | 1.260 | 1.240 | |

Worst | 1.987 | 1.289 | 1.314 | |

Average | 1.746 | 1.279 | 1.261 | |

Standard deviation | 0.162 | 0.010 | 0.024 | |

Coefficient of variation | 0.093 | 0.008 | 0.019 |

*x*as the reference point to identify the promising areas near the best solution. The parasitism operator introduces further merit to the computational process. Making changes in all the dimensions rather than in a small number of dimensions adds large perturbation to the ecosystem that maintains diversity while preventing premature convergence. Furthermore, due to the highly random nature of parasitism its uniquely produced solutions might be located in totally different areas. Hence, the operators' characteristics play pivotal roles in the successful performance of SOS. Figure 8 represents the volume of water released from the reservoir, the generated power, and the variation of reservoir storage during the operational period, respectively. Unlike the GA-driven results, the SOS algorithm's results were remarkably similar to the optimized values of NLP graphed in Figure 8. For instance, in Figure 8 during the period of the 30th to 42nd month, there was a low compatibility between the GA and the NLP results. It is seen in Figure 8 that the lowest compatibility between the GA and NLP occurred between periods 40 through 60 and between periods 1 through 10, respectively. As can be seen in Figure 8, the compatibility of WCA with NLP results is better than the GA's, yet it is not as efficient as the SOS.

_{best}## CONCLUDING REMARKS

This work introduced the SOS algorithm to the optimization of reservoir system operation and applied it to the Karun4 reservoir system with hydropower purpose, and to the optimal operation of a four-reservoir system. Our results from these two optimization problems indicate that the SOS algorithm outperformed the GA and a relatively new EA called the WCA. In the single reservoir problem with global optimal value equal to 1.213, the SOS algorithm, the WCA, and the GA determined 1.240, 1.260, and 1.535 as the optimal solutions, respectively. The superiority of the SOS algorithm over the other EAs was repeated in the four-reservoir optimizing problem with global optimal value equal to 308.29. In this regard, the GA, the WCA, and the SOS algorithm in their best performance among 10 runs converged to 97.46%, 99.56%, and 99.86% of the global optimal solution. Besides closeness of the best performance result to the global optimal solution calculated with the SOS algorithm, other practical merits are germane to the SOS algorithm. In comparison with the GA the SOS algorithm did not exhibit premature convergence, it had lower standard deviation, and had lower and more stable coefficient of variation. Moreover, unlike the WCA and the GA, the SOS does not require calibration of algorithmic parameters. The results showed that despite the simple principles of the SOS algorithm, it is capable of dealing with complexities and constraints associated with multi-reservoir operation optimization problems. This renders the SOS algorithm an attractive solution method of optimization problems in water resources management.

The superiority of the SOS algorithm over other EAs established in this study should not be considered applicable to all optimization problems. This paper's results are in line with the aforementioned ‘no free-lunch’ theorem, which emphasizes that a particular EA cannot optimally solve all well-posed optimizing problems. However, the SOS algorithm has unquestionable advantages, such as the simplicity of parameter specification, adding a perturbation to the search domain, and substituting a solution by evaluating the difference between other solutions that render it a very attractive EA. Future applications of the SOS algorithm could include solutions to conflicting objectives that appear in water resources management. For this purpose, multi-objective algorithm benchmarking can be accomplished with the NSGA-II, Borg, AMALGAM, and other algorithms which have shown good performance in solving multi-objective optimization problems.