Assessment of Water Resources Management Strategy Under Di ﬀ erent Evolutionary Optimization Techniques

: Competitive optimization techniques have been developed to address the complexity of integrated water resources management (IWRM) modelling; however, model adaptation due to changing environments is still a challenge. In this paper we employ multi-variable techniques to increase conﬁdence in model-driven decision-making scenarios. Here, water reservoir management was assessed using two evolutionary algorithm (EA) techniques, the epsilon-dominance-driven self-adaptive evolutionary algorithm ( ε -DSEA) and the Borg multi-objective evolutionary algorithm (MOEA). Many objective scenarios were evaluated to manage ﬂood risk, hydropower generation, water supply, and release sequences over three decades. Computationally, the ε -DSEA’s results are generally reliable, robust, e ﬀ ective and e ﬃ cient when compared directly with the Borg MOEA but both provide decision support model outputs of value.


Introduction
Water resource management problems (i.e., surface and groundwater) are complex due to their non-linear, dynamic, multimodal properties that need robust methods to solve, such as optimization algorithms [1] based on evolutionary algorithms (EAs) inspired from evolution and the natural selection of species [2,3]. Many of these have been proposed by researchers with different techniques, such as: the non-dominated sorting genetic algorithm (NSGA II) [4], multi-objective evolutionary algorithm based on decomposition (MOEA/D) [5], indicator-based evolutionary algorithm (IBEA) [6] and differential evolution (DE) [7]. Furthermore, approaches based on swarm intelligence include particle swarm optimization (PSO) [8] and ant colony optimization (ACO) [9] while the annealing process in metallurgy inspired simulated annealing (SA) [10]. A review of EAs and other metaheuristic algorithms and their applications can be found in [11,12]. Examples using these techniques for solving water resources management problems include Hurford et al., [13], and others [14][15][16][17] using ε-NSGA-II, MOEA/D, Borg MOEA and NSGA-II, respectively to optimize reservoir management strategy based on multidisciplinary objectives like flood control, hydropower generation, and water supply.
Benchmark functions such as DTLZ and WFG series were often used in comparative studies to assess algorithms' performance, as in [18][19][20][21], however they consider forward and easy to solve versus real-world problems [22]. These algorithms often have many parameters that require calibration, which has a major impact on computational performance and optimal achievement [23,24]. Karafotias et al., [25] presented a review of the approaches for the calibration and control of parameters. There are two types of EA parameter-setting problems categorised as (a) parameter tuning and (b) parameter 1.
In a minimization problem, a vector u = (u 1 , . . . , u M ) T is said to dominate another vector v = (v 1 , . . . , v M ) T if u i ≤ v i for i = 1, . . . , M and u v. This property may be denoted as u ≺v.

2.
A feasible solution x∈ X is called a Pareto-optimal solution, if there is no alternative solution y∈ X such that F(y) ≺ F(x).

3.
The Pareto-optimal set, PS, is the union of all Pareto-optimal solutions, and may be defined as PS = {x ∈ X : y ∈ X, F(y) ≺ F(x)}. 4.
The Pareto-optimal front, PF, is the set comprising the Pareto-optimal solutions in the objective space. It may be expressed as PF = {F(x)|x ∈ PS}.

Details of Epsilon-Dominance-Driven Self-Adaptive Evolutionary Algorithm (ε-DSEA) Optimization Algorithm
The algorithm is based on the main principles of multi-objective evolutionary algorithms (MOEAs) e.g. recombination, mutation and dominance sorting. However, novel techniques are included to enhance the algorithm's ability to handle the complexities of different problem environments. These techniques are: 1.
Diversity expansion to increase decision variables' search space exploitation 2.
Self-adaptive operators' parameters for parameters in process tuning 3.
Exploration extension for algorithm revival and stagnation coping 4.
Virtual dominance archive to improve diversity and convergence.

Diversity Expansion
The search procedure in an optimization algorithm has two main components, exploration and exploitation. Evidence in the literature indicates the best results are achieved if exploration and exploitation are deployed preferentially in the early and latter stages of the search, respectively [47,48]. Accordingly, a procedure that safeguards diversity in the population at the start is incorporated in the proposed algorithm which employs all the available recombination operators at the initial stage. After the initial random seeding, the algorithm uses each recombination operator to generate new offspring, selecting parents from the entire population. If more parents are needed (e.g., in case of odd number of parents), they are selected from the population using a binary tournament selection. Figure 1 illustrates the procedure by which the parents are selected. Illustrates operator's parents' selection from the entire population candidates after the initial random seeding at the beginning of the evaluation process.

Self-Adaptive Mechanism and Formulae
In each generation, the recombination operators are selected on a competitive basis, according to the proportion of dominance solutions in the archive (NDS) contributed by each operator. Thus, the selection probabilities for the recombination operators are obtained as follows [30,49].
where is the probability of i th recombination operator, NDSi is the number of solutions in the archive contributed by the i th recombination operator, NRO is the number of recombination operators; The constant 1.0 is used to avoid probability values of zero.
However, operator's dominance achievement is sensitive to the relevant parameter setting. This problematic, (e.g., parameter control problem) was classified into three categories depending on the way the parameter variation is accomplished [26]: (a) deterministic, (b) adaptive, and (c) selfadaptive. Deterministic control is based on rules that are specified a priori [50,51]. In self-adaptive control, the parameters may be encoded to evolve in the genotype such that, for example, mutation and recombination are applied to the decision variables also [52,53], which extends the search space to cover the parameter values and consumes more time during the optimization processes [26]. The adaptive method (b) was considered more effective in solving complex problems, as a feedback from the optimization process is set to adjust parameter values during optimization progresses [26,54]. The aforementioned technique was adopted by many researchers to improve algorithms' performance, as in [30,49,55−57], however none of these (and others) develop a self-adaptive technique that is sensitive to optimality achievement during evaluation progress [31].
The success of any operator depends on the chosen values of the parameters that directly affect its performance. Any operator may lead an algorithm to suboptimal solutions because of unsatisfactory parameter calibration. However, parameter calibration is extremely challenging. This difficulty provided the motivation for establishing a dynamic relationship between the values of the control parameters of the recombination operators and their relative effectiveness, to obviate the need for fine tuning. The efficiency of the optimization algorithm is thus improved by continuously seeking to improve the collective effectiveness of the recombination operators. In other words, the formulation developed herein allows the values of the control parameters of each recombination operator to improve adaptively based on the success of the recombination operator compared to the rest of the recombination operators. Illustrates operator's parents' selection from the entire population candidates after the initial random seeding at the beginning of the evaluation process.

Self-Adaptive Mechanism and Formulae
In each generation, the recombination operators are selected on a competitive basis, according to the proportion of dominance solutions in the archive (NDS) contributed by each operator. Thus, the selection probabilities for the recombination operators are obtained as follows [30,49].
where P NDS i is the probability of i th recombination operator, NDS i is the number of solutions in the archive contributed by the i th recombination operator, NRO is the number of recombination operators; The constant 1.0 is used to avoid probability values of zero.
However, operator's dominance achievement is sensitive to the relevant parameter setting. This problematic, (e.g., parameter control problem) was classified into three categories depending on the way the parameter variation is accomplished [26]: (a) deterministic, (b) adaptive, and (c) self-adaptive. Deterministic control is based on rules that are specified a priori [50,51]. In self-adaptive control, the parameters may be encoded to evolve in the genotype such that, for example, mutation and recombination are applied to the decision variables also [52,53], which extends the search space to cover the parameter values and consumes more time during the optimization processes [26]. The adaptive method (b) was considered more effective in solving complex problems, as a feedback from the optimization process is set to adjust parameter values during optimization progresses [26,54]. The aforementioned technique was adopted by many researchers to improve algorithms' performance, as in [30,49,[55][56][57], however none of these (and others) develop a self-adaptive technique that is sensitive to optimality achievement during evaluation progress [31].
The success of any operator depends on the chosen values of the parameters that directly affect its performance. Any operator may lead an algorithm to suboptimal solutions because of unsatisfactory parameter calibration. However, parameter calibration is extremely challenging. This difficulty provided the motivation for establishing a dynamic relationship between the values of the control parameters of the recombination operators and their relative effectiveness, to obviate the need for fine tuning. The efficiency of the optimization algorithm is thus improved by continuously seeking to improve the collective effectiveness of the recombination operators. In other words, the formulation developed herein allows the values of the control parameters of each recombination operator to improve adaptively based on the success of the recombination operator compared to the rest of the recombination operators. Table 1 shows the lower and upper bounds of the operator control parameters. If an operator's ability to contribute offspring to the dominance archive is decreased, its selection probability P NDS i will decrease according to Equation (2). In turn, the values of the relevant control parameter decrease and the recombination operator's ability to contribute new offspring to the archive will improve. It is worth noting that, initially, all the recombination operators have an equal selection probability (P NDS i ) of 1/NRO. During the evaluation process the P NDS i value for any recombination operator changed along with its control parameters, according to its contribution in the dominance archive. If any recombination operator is relatively unsuccessful, its selection probability (P NDS i ) and parameters controls will decrease. If the effectiveness of another recombination operator decreases, the selection probabilities of some or all the other recombination operators will increase with the values of their control parameters. In this way a dynamic equilibrium is maintained among the operators' selection probabilities, which in turn regulates the operator control parameters. A set of formulae (equations) developed ensemble parameters' tuning with the relevant dominance attainment. The parameters' tuning domains (i.e., tuning range) were set based on default or recommended values suggested in the literature, and experimental investigation carried out on common test functions, as illustrated in Table 1. Figure 2a illustrates the relationships between operators' dominance attainment and their control parameters values and shows how operators' parameters auto-tuned according to the operator successful to produce non-dominated solutions in the dominance archive.

Exploration Extension Mechanism
This mechanism is based on initializing (resetting) all operators' selection probabilities P NDS i uniformly to1/NRO. It aims to provide an equal opportunity for all the operators, by assessing the performance best on the most recent results. Otherwise, the previously successful operators with more solutions in the archive would continue to dominate based on past performance as dictated by Equation (2).
The number of resets depends on a random integer N r such that N r ∈ N + ∈ [1,3]. When the algorithm starts, an N r value is selected at random and the maximum permissible number of function evaluations NFE max is divided by N r + 1 to determine the reset interval E r . For example, if NFE max = 300,000 and N r = 2, the reset occurs at every E r = 100,000 function evaluations. Hence, in this case, two resets occur during the entire optimization. Formally, where E r is the reset interval. Figure 2b shows an example of the resetting process and its relation with self-adaptive mechanism to extend algorithm explorations and escaping from possible local optima.
Water 2019, 11, x FOR PEER REVIEW 6 of 23 Figure 2b shows an example of the resetting process and its relation with self-adaptive mechanism to extend algorithm explorations and escaping from possible local optima.

Virtual Dominance Archive
In early stages of an evaluation process for constraints problems with enormous decision variables, the ε-dominance archive techniques (Section 2.1.3) tend to maintain only the nondominated solutions in the dominance archive. Experimental tests on such problems show only one non-dominated solution maintained in the archive while exploring the design space for feasible solutions. Hence, the operators' parameters will be on its minimum values during this stage in the evaluation process using the proposed self-adaptive mechanism. To overcome this issue, a virtual

Virtual Dominance Archive
In early stages of an evaluation process for constraints problems with enormous decision variables, the ε-dominance archive techniques (Section 2.2.3) tend to maintain only the non-dominated solutions in the dominance archive. Experimental tests on such problems show only one non-dominated solution maintained in the archive while exploring the design space for feasible solutions. Hence, the operators' parameters will be on its minimum values during this stage in the evaluation process using the proposed self-adaptive mechanism. To overcome this issue, a virtual dominance archive was developed by randomly generating a virtual number of the dominance solutions for the selected operator to preserve diversity and early convergence exploration for feasible solutions using the entire parameter's domain.

Constraint Handling Strategy
The fundamentals of evolutionary algorithms are based on handling only unconstraint optimization problems [35], many techniques were proposed for constraint problems, like penalty function, special representations and operators, and repair method [58]. Here, the penalty function technique is adopted as follows [59]: where F (x) is the expanded objective function, and P(x) is the constraint violation amount, which can be expressed as (based on Equation (1)): where A i and B j ∈ R + , are penalty factors. I and J are the total numbers of inequity and equity constraints, respectively.
MOEAs' effectiveness is commonly measured using quantitative metrics like the hypervolume metric [65] which evaluate the non-dominated solutions' hypervolume, and generational distance metric [66] which measure the average distance between the dominance solutions and the closer Pareto-front set. However, these metrics (and others) may provide misguiding results and most of their design principles depends on the true Pareto-front, which is unknown in real-world problems [22].
Accordingly, the comparative assessment of ε-DSEA are based on a real-world engineering problem. Here, the state-of-the-art Borg MOEA [30] was adopted for comparative purposes since it outperforms or met other state-of-the-art algorithms' achievement, such as: ε-MOEA, ε-NSGA-II, MOEA/D, GBE3, OMOPSO, IBEA, NSGA-II, AMALGAM [18,30,[67][68][69][70]. Borg employs many MOEAs' design principles based on previous works like; recombination, mutation, and dominance sorting (e.g., ε-box). The authors present novel techniques to improve the exploration and exploitation process including; ε-progress indicator of stagnation and improvement, population expansion to preserve diversity exploration, multiple recombination operators for search variations, and self-adaptive of operator. A concise detail of these techniques are presented below, more details are presented in aforementioned literatures.
Borg MOEA uses an active population of solutions and an external archive that stores dominant solutions, and the population size is proportional to the archive size. Initially, the archive is empty; hence an initial population size is required. Subsequently, the population size changes as follows [30] where NP and NA are the population and archive sizes, respectively, and γ is the ratio of the population size to the archive size, and equal to 4 [30]. The ε-progress index measures the improvements while searching for new solutions. If the algorithm finds new dominant solutions in a new unoccupied ε-box (if the new dominant solutions have different ε-box indices) it means there is improvement, otherwise no improvement flag will mark a stagnation sign. If the last case continues for a number of evaluations, a revival process named "restart" will be triggered to escape from possible local optima. The restart involves emptying the population and re-populating based on the population to archive ratio (Equation (6)). The population is refilled using all solutions in the archive. Any remaining empty slots in the population are filled with solutions created by uniform mutation of solutions that are selected randomly from the archive.
The trigger for the revival process depends on any of the following three conditions: a.
If there is no change in the archive size for a certain number of evaluations; b.
If there is no improvement indicated by the -progress indicator; and c.
If the current population to archive ratio exceeds 1.25×γ Borg follows the same crossover and mutation techniques mentioned in Section 2.2, and employs Equation (2) to self-adaptive operator's selection. In Borg, the relevant operators' parameters have fixed pre-execution values during evaluation process.

Identification of a Real-World Experimental Test Problem
A case study in Iraq's Diyala river basin was adopted as a real-world IWRM problem which is more complex than common benchmark test functions [22]. GWP, [71]. The Global Water Partnership defines the IWRM as "IWRM is a process which promotes the co-ordinated development and management of water, land and related resources, in order to maximize the resultant economic and social welfare in an equitable manner without compromising the sustainability of vital ecosystems". Authors and institutes adopt different water management concepts (about 41 variant possible explanations for the term "integrated") due to the generalization in IWRM definition. Some examples are: water supply and water demands; surface water and groundwater; water quantity and water quality; urban and rural water issues; government and NGOs (non-governmental organizations) [72,73]. The river basin has two multipurpose dams, Derbendikhan just at the northern international border in Sulaymaniya governorate, and Himren in the middle part of the basin in the Diyala governorate ( Figure 3). Here, Derbendikhan dam's operation strategy for the next three decades was selected as a benchmark problem. Based on monthly dataset from 1981 to 2012 (33 years), a total of 396 decision variables (reservoir releases) need to be managed during the time-scale. Generating hydropower is the main current operation target, hence power penstocks (tunnels) are the main reservoir outlet of the proposed management model. governorate, and Himren in the middle part of the basin in the Diyala governorate ( Figure 3). Here, Derbendikhan dam's operation strategy for the next three decades was selected as a benchmark problem. Based on monthly dataset from 1981 to 2012 (33 years), a total of 396 decision variables (reservoir releases) need to be managed during the time-scale. Generating hydropower is the main current operation target, hence power penstocks (tunnels) are the main reservoir outlet of the proposed management model.

Objectives Functions Formulae
The reservoir water budget is governing by the water balance equation, as:

Objectives Functions Formulae
The reservoir water budget is governing by the water balance equation, as: where S D t and S D t+1 are the reservoir storage at time t and t+1, I D t and R D t are reservoir inflows and releases, respectively. E D t is the evaporation losses from reservoir surface, P D t is the direct rainfall on the reservoir. While, SE D t and GR D t are seepage losses and groundwater recharges from the reservoir, respectively.
The reservoir operation strategy (F D ) is represented by the following multi-objective (or many for more than 3 objectives) formula: where f winterD is for maximizing winter storage to fulfil summer demands, f summerD is for minimizing summer storage to absorb expected flood wave in the next season, f powerD is for maximizing hydropower generation, f Del−SW is for minimizing agriculture projects' water deficit, and f regD is for minimization releases fluctuation. These targets represent the following aspects: social ( f winterD and f summerD ); economic ( f powerD , f Del−SW ); and environmental ( f regD ). The details of these objectives functions are as follows: Where: S D max = maximum allowable reservoir storage S D minp = minimum allowable reservoir storage for hydropower generation T W , T S and T = winter, summer and total operation periods, respectively. Pw D t = hydropower generation at time t Pw D max = maximum hydropower generation PD t = projects' water demands at time t PD max = maximum projects' water demands Del M t = delivered water at time t R D max = maximum reservoir releases at time t C P = penalty factor includes all the violations of the model, which could be expressed NC = number of constraints g i = penalty function for the (i th ) constraint A = a positive real number The hydropower generation can be expressed as: where (Q tuD t ) is the turbine discharge, (H nD t ) is the net head between reservoir level and the tail water after the power plant, (η D e ) is the efficiency of power plant, and (γ w ) is the water density.

Reservoir System Constraints
The reservoir storage is limited between the minimum and maximum allowable storage, 283.48 ≤ S D Max 0, 249000 − Pw D t (20)

Computational Properties
The computational parameters of the problems were 2.0 × 10 6 function evaluations and ε = 0.1 for three objectives, and ε = 0.5 for five objectives, with 10 and 20 runs for both algorithms. The minimum population size was 100 while the maximum was 1000. A Dell OptiPlex 780 computer was used (Core Duo 2 E8400, 2 × 3.0 GHz, 8.0 GB RAM, Ubuntu 16.04 operating system). Table 2 shows the parameter values used for both algorithms. A program (code) in C language was developed to build the current model. L is the number of decision variables. The permissible range for dynamic parameters is shown in brackets. The parameters σ η and σ ζ are defined in Table 1. a The initial values of dynamic parameters used in ε-DSEA are as shown for Borg MOEA. Figure 4 illustrates the Pareto-front for 20 replicated random runs of the case study benchmark problem for both algorithms using three objectives. Although both algorithms converged to possible optima, ε-DSEA shows better reliability. Notably, Borg MOEA has faced some challenges in six trials, as it was stagnant in four and had overdue convergence in two. Some of these contain only one solution, as highlighted with dotted lines. In addition, some others have discontinuous Pareto-front (with gaps) in the objective search space. This behaviour reduces an algorithm's reliability to produce near optimal solutions over replicated random execution (for example when testing the confidence of the model output), and is a key factor when solving more complex problem using high-performance computer resources (e.g., parallel processing with multi-core). Conversely, ε-DSEA shows reliability over 20 runs to converge to the possible optimal solution, with commonly continuous Pareto-front. Hence, less randomness creeps into runs increasing the confidence that optimal solutions can be achieved.

Algorithms' Reliability
(with gaps) in the objective search space. This behaviour reduces an algorithm's reliability to produce near optimal solutions over replicated random execution (for example when testing the confidence of the model output), and is a key factor when solving more complex problem using highperformance computer resources (e.g., parallel processing with multi-core). Conversely, ε-DSEA shows reliability over 20 runs to converge to the possible optimal solution, with commonly continuous Pareto-front. Hence, less randomness creeps into runs increasing the confidence that optimal solutions can be achieved.  A high-dimension problem (5 objectives) was employed for advance algorithms' assessment as shown in Figure 5. The optimum solutions' median for 10 trials of both algorithms is presented. Notably, both algorithms produce possible optimum solutions, but the ε-DSEA has slightly more reliable trends over execution repetitions as supported by the self-adaptive parameters' technique used in this EA, which was also approved by [74] for 20 runs. Insight investigation shows that a Borg MOEA stack in local optima twice (run 3 and 7), and adapt with PCX operator for all trials, as in three objectives scenario; while ε-DSEA adapts at the initial stage with SBX operator, then with PCX and SPX operators in parallel for the rest of evaluation process ( Figure S1 in the supplementary data). The resetting technique's effectiveness is obvious over changing the trend of the operators' adaptation to escape from a local optima pitfall. Execution trial No. 4 shows competitive achievement of both algorithm, which may consider for comparative investigation. A high-dimension problem (5 objectives) was employed for advance algorithms' assessment as shown in Figure 5. The optimum solutions' median for 10 trials of both algorithms is presented. Notably, both algorithms produce possible optimum solutions, but the ε-DSEA has slightly more reliable trends over execution repetitions as supported by the self-adaptive parameters' technique used in this EA, which was also approved by [74] for 20 runs. Insight investigation shows that a Borg MOEA stack in local optima twice (run 3 and 7), and adapt with PCX operator for all trials, as in three objectives scenario; while ε-DSEA adapts at the initial stage with SBX operator, then with PCX and SPX operators in parallel for the rest of evaluation process ( Figure S1 in the supplementary data). The resetting technique's effectiveness is obvious over changing the trend of the operators' adaptation to escape from a local optima pitfall. Execution trial No. 4 shows competitive achievement of both algorithm, which may consider for comparative investigation. Hence, the proposed mechanism provides advance diversity and balancing between exploration and exploitation process toward possible Pareto-front set.

Algorithms' Robustness and Efficiency
In ε-DSEA, non-dominated feedback loops control the operators' adaptation and their parameters. Figure 6a illustrates the self-adaptive operators' parameter-tuning behaviour of ε-DSEA during the evaluation process. The most effective operators adopted to generate dominance solutions for the best trial are SBX, PCX, and SPX. Initially the virtual dominance archive mechanism tuned the Hence, the proposed mechanism provides advance diversity and balancing between exploration and exploitation process toward possible Pareto-front set.

Algorithms' Robustness and Efficiency
In ε-DSEA, non-dominated feedback loops control the operators' adaptation and their parameters. Figure 6a illustrates the self-adaptive operators' parameter-tuning behaviour of ε-DSEA during the evaluation process. The most effective operators adopted to generate dominance solutions for the best trial are SBX, PCX, and SPX. Initially the virtual dominance archive mechanism tuned the operator's parameters when only one solution is kept in the dominance archive. Then the SBX operator was adopted until the first resetting trigger at 5.0 × 10 5 function evaluation. The PCX operator is then involved by increasing the variation parameters (σ η and σ ζ ) to about 0.15. The SPX operator is also involved at the same time when its parameter (λ) changed to about 2.7. Both PCX and SPX operators compete to explore dominance solutions till the third resetting trigger, and after that the SPX operator starts to generate more dominance solutions in the dominance archive. Increasing PCX and SPX parameters will generate new offspring farther away from their original parents, which will increase algorithm exploration in the design search space. The algorithms' convergence (efficiency) also investigated using the decision variables vector ( ) development in the dominance archive during the evaluation process. The is equal to + + + ⋯ + , where x1 to xn are the decision variables. Based on the best solution achieved, Figure 6b shows convergence of both algorithms. Both achieved early convergence, but ε-DSEA converged faster, hence ε-DSEA's efficiency was endorsed in the proposed test problem.
The progress of the objectives' convergence of both algorithms over 10 iterations of highdimension problem is presented in Figure S2 and S3 in the supplementary data, respectively. Early convergence was achieved by both algorithms, ε-DSEA converged at 1.25 × 10 4 function evaluations for all iterations, and Borg MOEA converged at 25 × 10 4 . The ε-DSEA needs less execution time to achieve solutions. Where there are limited computational resources (e.g., CPU, Ram, etc.,) this achievement is significant. Furthermore, Borg MOEA suffered significant and interim stagnation in 7 trials (2, 3, 7, 9, and 4, 6, 10, respectively) in the early stage of evaluations. Only three out of 10 trials maintained dominance solutions improvement over the entire evaluation. The PCX operator's adaptation with fixed parameters and recycling repetitively archive's dominance solutions may restrict the extent of the algorithm's exploration in the design search space. Conversely, only one trial (no. 9) suffered significate stagnation in ε-DSEA, however the expansion diversity and resetting techniques succeed in reviving the algorithm's exploration to find new dominance solutions in the dominance archive. The robustness of ε-DSEA to escape from local optima are evident. Figure 7 shows trial no. 4 as a sample of convergence progress, since both algorithms achieved competitive solutions (based on Figure 5). The algorithms' convergence (efficiency) also investigated using the decision variables vector (X dv ) development in the dominance archive during the evaluation process. The X dv is equal to , where x 1 to x n are the decision variables. Based on the best solution achieved, Figure 6b shows X dv convergence of both algorithms. Both achieved early convergence, but ε-DSEA converged faster, hence ε-DSEA's efficiency was endorsed in the proposed test problem.
The progress of the objectives' convergence of both algorithms over 10 iterations of high-dimension problem is presented in Figures S2 and S3 in the supplementary data, respectively. Early convergence was achieved by both algorithms, ε-DSEA converged at 1.25 × 10 4 function evaluations for all iterations, and Borg MOEA converged at 25 × 10 4 . The ε-DSEA needs less execution time to achieve solutions. Where there are limited computational resources (e.g., CPU, Ram, etc.,) this achievement is significant. Furthermore, Borg MOEA suffered significant and interim stagnation in 7 trials (2, 3, 7, 9, and 4, 6, 10, respectively) in the early stage of evaluations. Only three out of 10 trials maintained dominance solutions improvement over the entire evaluation. The PCX operator's adaptation with fixed parameters and recycling repetitively archive's dominance solutions may restrict the extent of the algorithm's exploration in the design search space. Conversely, only one trial (no. 9) suffered significate stagnation in ε-DSEA, however the expansion diversity and resetting techniques succeed in reviving the algorithm's exploration to find new dominance solutions in the dominance archive. The robustness of ε-DSEA to escape from local optima are evident. Figure 7 shows trial no. 4 as a sample of convergence progress, since both algorithms achieved competitive solutions (based on Figure 5). 7 trials (2, 3, 7, 9, and 4, 6, 10, respectively) in the early stage of evaluations. Only three out of 10 trials maintained dominance solutions improvement over the entire evaluation. The PCX operator's adaptation with fixed parameters and recycling repetitively archive's dominance solutions may restrict the extent of the algorithm's exploration in the design search space. Conversely, only one trial (no. 9) suffered significate stagnation in ε-DSEA, however the expansion diversity and resetting techniques succeed in reviving the algorithm's exploration to find new dominance solutions in the dominance archive. The robustness of ε-DSEA to escape from local optima are evident. Figure 7 shows trial no. 4 as a sample of convergence progress, since both algorithms achieved competitive solutions (based on Figure 5).

Algorithms' Effectiveness
For real-world multi-objective problems, and especially in water resources management problems, the true Pareto-front (e.g., optimum solution) is unknown [22], and it is difficult to measure an algorithm's effectiveness for such problems, as other relevant factors should also be evaluated such as the coverage of the Pareto-front and its extent in the objective space [34]. Hence, the qualitative comparison was often based on the best solution achieved over several replicated trials (e.g., equal or more than 20 runs). The results here show the reliability of both EA models but better computational performance by ε-DSEA. Table 3 demonstrates results' analysis of both algorithms' achievement based on best optimum solution to maximise hydropower generation, as it is one of the main dam's operation targets. The gross sum of hydropower, storage, and releases of the reservoir were presented to demonstrate the contrast between two algorithms' achievement, based on the relevant optimization techniques. The results of both are harmonic, with advance merit of ε-DSEA. Figure 8a,b depict algorithms' attainment of 10 multi-objective multi-variable trials and show the consistency of the ε-DSEA is marginally better than Borg notably from the period 1 to 216 months. The same behaviour also achieved for the next period, which reflect ε-DSEA algorithms' ability to generate possible competitive optimal solution with fewer replicated trials.

Strategic Achievement
Consistency with the relevant historical (actual) dataset should also be reviewed during decision making trials (i.e. solutions' quality), as in Figure 8c,d. Both algorithms achieved competitive results, but ε-DSEA's result has better reaction to flood waves and better agreement with the historical data. Here, spillway discharge did not factor in the models or the developed management model, since flood waves events usually last hours or a couple of days in the investigated region. Accordingly, no relevant sensitive reaction was observed since monthly average management was adopted by the model.  Figure 8a,b depict algorithms' attainment of 10 multi-objective multi-variable trials and show the consistency of the ε-DSEA is marginally better than Borg notably from the period 1 to 216 months. The same behaviour also achieved for the next period, which reflect ε-DSEA algorithms' ability to generate possible competitive optimal solution with fewer replicated trials.
Consistency with the relevant historical (actual) dataset should also be reviewed during decision making trials (i.e. solutions' quality), as in Figure 8 c,d. Both algorithms achieved competitive results, but ε-DSEA's result has better reaction to flood waves and better agreement with the historical data. Here, spillway discharge did not factor in the models or the developed management model, since flood waves events usually last hours or a couple of days in the investigated region. Accordingly, no relevant sensitive reaction was observed since monthly average management was adopted by the model. As a source of renewable energy, hydropower generation is one of the key-operational targets of the tested real-world problem, and often for any multipurpose dam projects, that needs to be carefully management under different operation scenarios, such as flood risk management. Two objectives were adopted, and , the later selected for insight investigation, as it is the most critical operation scenario that may affect other operational targets. Figure 9a,b illustrate solution distribution density achieved by both algorithms over the power generation domain. In general, ε-DSEA achieved high repetition of 30 to 50 Mw solutions, and gentle gradient repetition after that, while Borg MOEA has steeper gradient repetition starting from 30 Mw and thereafter. ε-DSEA achievement offers insight for investment decision making as minimum power generation of 30 Mw could be guaranteed for next three decades.
Hydropower generation depends on two variables, turbine's discharge and water net head, as in Equation 15 (turbine's efficiency and water's specific weight assumed constant). Figure 9c,d illustrate hydropower generation solution's space achieved by both algorithms. Competitive As a source of renewable energy, hydropower generation is one of the key-operational targets of the tested real-world problem, and often for any multipurpose dam projects, that needs to be carefully management under different operation scenarios, such as flood risk management. Two objectives were adopted, f winterD and f summerD , the later selected for insight investigation, as it is the most critical operation scenario that may affect other operational targets. Figure 9a,b illustrate solution distribution density achieved by both algorithms over the power generation domain. In general, ε-DSEA achieved high repetition of 30 to 50 Mw solutions, and gentle gradient repetition after that, while Borg MOEA has steeper gradient repetition starting from 30 Mw and thereafter. ε-DSEA achievement offers insight for investment decision making as minimum power generation of 30 Mw could be guaranteed for next three decades.
Water 2019, 11, x FOR PEER REVIEW 2 of 23 tends to maximize power generation by releasing more water, but solutions are irregularly deployed (clustering) which possibly due to local optima pitfall. For example, in the region of <80m net head ( ) and >600 MCM releases (Figure 9c,d), only five solutions formed as two groups achieved (green colour), which is not the case in ε-DSEA.
Recreation is another target to optimize for revenue. Figure 10 shows reservoir surface area achieved by both algorithms for the considered time-scale of scenario. The mean and median surface area were about 52 and 49 km 2 , 51 and 44 km 2 achieved by ε-DSEA and Borg respectively. The small violation between these values indicates competitive results' distribution, corollary more solutions greater than these values achieved by ε-DSEA. Hence, projects' revenues could be improved even in such a critical scenario. Hydropower generation depends on two variables, turbine's discharge and water net head, as in Equation (15) (turbine's efficiency and water's specific weight assumed constant). Figure 9c,d illustrate hydropower generation solution's space achieved by both algorithms. Competitive distribution over solution space was accomplished by both algorithms for turbine discharge ≤600 MCM, while better space's exploitation for greater values (>600 MCM) achieved by ε-DSEA. Borg tends to maximize power generation by releasing more water, but solutions are irregularly deployed (clustering) which possibly due to local optima pitfall.
For example, in the region of <80m net head (H nD ) and >600 MCM releases (Figure 9c,d), only five solutions formed as two groups achieved (green colour), which is not the case in ε-DSEA.
Recreation is another target to optimize for revenue. Figure 10 shows reservoir surface area achieved by both algorithms for the considered time-scale of f summerD scenario. The mean and median surface area were about 52 and 49 km 2 , 51 and 44 km 2 achieved by ε-DSEA and Borg respectively. The small violation between these values indicates competitive results' distribution, corollary more solutions greater than these values achieved by ε-DSEA. Hence, projects' revenues could be improved even in such a critical scenario.
five solutions formed as two groups achieved (green colour), which is not the case in ε-DSEA.
Recreation is another target to optimize for revenue. Figure 10 shows reservoir surface area achieved by both algorithms for the considered time-scale of scenario. The mean and median surface area were about 52 and 49 km 2 , 51 and 44 km 2 achieved by ε-DSEA and Borg respectively. The small violation between these values indicates competitive results' distribution, corollary more solutions greater than these values achieved by ε-DSEA. Hence, projects' revenues could be improved even in such a critical scenario.

Algorithms' Optimization Techniques
The technique of Borg MOEA tends to adapt based on operator experience after finding possible feasible solutions. The SBX operator was adopted in the early stage of evaluation process, then PCX operator adopted to the end. Zheng et al., [48] observed that, for two-objective problems, Borg MOEA tended to converge prematurely and population diversity decreased relatively rapidly. This is because Borg MOEA does not maintain a separate transient sub-population of offspring as in NSGA-II. Offspring that dominate its parents immediately replaces one of the parents; the choice of the parent that is replaced is random. As new solutions are introduced, the selection pressure on less competitive solutions increases, due to the binary tournament selection used for a crossover. Fitting solutions have a higher probability of selection for a crossover, leading to more exploitation and less exploration and thus less diversity. Secondly, the injection trigger that depends mainly on ε-progress indicators, did not always succeed in reflecting stagnation during the evaluation process. Thirdly, PCX operator produces offspring in the vicinity of the parents. If the PCX operator creates solutions around the best solutions, the PCX-generated solutions quickly dominate the archive, leading to more exploitation, less exploration and consequently relatively rapid loss of diversity. As stated previously, the recombination operators are deployed in proportion to the number of offspring they contributed in the archive.
In Borg MOEA the operator that produces more successful (i.e., non-dominated) offspring is deployed more frequently. However, as the search progresses and the balance between exploration and exploitation shifts gradually, it is desirable that the operators be deployed based on the current status of the search rather than their previous performances or cumulative successes. In other words, the selection of the operators should recognize the current performance also. Hence, the proposed performance assessment of the operators relies on the results from the current phase of the search rather than the cumulative performance to date.
The novel technique of virtual dominance archive used in ε-DSEA removes operator bias and extends the algorithm's exploration by tuning the operators' parameters control within the specified domain, resulting in a robust convergence progress at initial stage of the evaluation process. In the same context, resetting these parameters' values during evaluation process help the algorithm to escape from local optima pitfall, and improve its exploitation to find new non-dominated solutions.

Water Resources Management Case Study
Although both algorithms achieved possible optimum solutions, ε-DSEA generates more robust results. Based on gross sum of five objective problems (Table 3), an extra 520 MW and 46.31 MCM of hydropower generation and reservoir storage was achieved, respectively. The mean and median water head achieved by ε-DSEA were about 1.5 m higher than those of Borg MOEA in all cases. This will provide advance security against possible dam failures in the future as the water head should be above 455.0 m.a.s.l [75]. Although this problematic was not considered in the current model as an objective or a constraint, all the relevant mean and median values achieved by both algorithms satisfy this restriction. This area of study should be considered in future work. Figure 11a,b illustrates results' quality of reservoir water level (m.a.s.l) produced by both algorithms over 10 runs, based on maximizing hydropower generation. Notably, more results of ≥460.0 and fewer of ≤440.0 were generated by ε-DSEA. The same achievement observed in the critical scenario of minimizing storage in summer, represented by f summerD , as in Figure 11c,d. Nevertheless, the latter has better results' than Borg MOEA regarding the current potential risk. Moreover, the harmony of results' distribution over 10 trials is evident, endorsing the reliability of ε-DSEA. the master decision variables' search space is relevant to reservoir releases while other dependent variables are calculated accordingly (e.g., power generation, storage, area, etc.,), the compete optimality achievement mapped to those sub-variables. This is not the case in other optimization algorithms, since often competitive investigation only covers the master decision variables and/or objectives' search space, either in test functions or in real-world case study benchmarks. Insight or deep diagnosis assessment considering indirect variables should be used in such competitive studies. The results suggest water resources decision makers are advised to implement different optimization algorithms, especially when solving multi or many objective problems, to explore possible new optimal results with advance quality, and to reinforce results' reliability. The ε-DSEA achievement was previously assessed using more complex water resources problem, as in [74,76], however, more assessment is recommended to solve different problem environments.

Conclusions
In this research, strategic planning of water resources under competitive optimization techniques was investigated. Decision makers often adopt optimization techniques to evaluate a wide range of competitive water resources management decisions. However, past studies demonstrate that an algorithm's optimal output varies according to the problem environment. Furthermore, slave (dependent) variables' quality are often not analysed in depth, which could play an important role in improving a project's economic success. Here, a comparative assessment of two optimization techniques was tested against a real-world water resources strategy. The Borg MOEA and the ε-DSEA performance was contrasted based on the relevant strategic plans using objective functions (i.e., the Pareto-front), master and slave variables.
The results by both models showed possible optima, with an advanced reliability and robustness when using ε-DSEA as it provided consistent results closer to near-optimal solutions. Both algorithms have auto-adaptive operator techniques, but Borg appoints a mono-operator after a specific number of evaluations (e.g. PCX), while multi-operators sequence during the evaluation stages (e.g. starting by SBX operator and ending with PCX and SPX operators in parallel). The Borg MOEA reviver's techniques should consider this drawback, as early stages stagnation is evident, especially in many- Within this framework, good exploration-exploitation balance of slave (or secondary) solution search space of hydropower generation and reservoir surface area were attained (Figure 9). Although the master decision variables' search space is relevant to reservoir releases while other dependent variables are calculated accordingly (e.g., power generation, storage, area, etc.), the compete optimality achievement mapped to those sub-variables. This is not the case in other optimization algorithms, since often competitive investigation only covers the master decision variables and/or objectives' search space, either in test functions or in real-world case study benchmarks. Insight or deep diagnosis assessment considering indirect variables should be used in such competitive studies.
The results suggest water resources decision makers are advised to implement different optimization algorithms, especially when solving multi or many objective problems, to explore possible new optimal results with advance quality, and to reinforce results' reliability. The ε-DSEA achievement was previously assessed using more complex water resources problem, as in [74,76], however, more assessment is recommended to solve different problem environments.

Conclusions
In this research, strategic planning of water resources under competitive optimization techniques was investigated. Decision makers often adopt optimization techniques to evaluate a wide range of competitive water resources management decisions. However, past studies demonstrate that an algorithm's optimal output varies according to the problem environment. Furthermore, slave (dependent) variables' quality are often not analysed in depth, which could play an important role in improving a project's economic success. Here, a comparative assessment of two optimization techniques was tested against a real-world water resources strategy. The Borg MOEA and the ε-DSEA performance was contrasted based on the relevant strategic plans using objective functions (i.e., the Pareto-front), master and slave variables.
The results by both models showed possible optima, with an advanced reliability and robustness when using ε-DSEA as it provided consistent results closer to near-optimal solutions. Both algorithms have auto-adaptive operator techniques, but Borg appoints a mono-operator after a specific number of evaluations (e.g., PCX), while multi-operators sequence during the evaluation stages (e.g., starting by SBX operator and ending with PCX and SPX operators in parallel). The Borg MOEA reviver's techniques should consider this drawback, as early stages stagnation is evident, especially in many-objectives problem. The ε-DSEA escapes from local optima by employing by-stage operators' parameters, which may be useful when applying to real-world problem solving.
The compete achievement of ε-DSEA was mapped onto the relevant water resources strategic plan. In all adopted scenarios for the real-world case study results show that extra hydropower, reservoir storage and surface area can be achieved. The releases have better consistency and sensitivity with the historical dataset and flood waves. The model outputs can be used to manage power generation to support, for example, an investment opportunity while still promising recreation investment opportunities achieved by maintaining larger reservoir surface area over the adopted time-scale. The results demonstrate the importance of insight and in-depth analysis of relevant objectives and variables using EA models.
The ε-DSEA and the relevant approach could be evolve for similar and/or even more complex real-world problems, such as groundwater management, water supply system, water allocation, etc., by adding and/or modifying objective functions, decision variables, and constraints.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/11/10/2021/s1: Figure S1: Active ε-DSEA operators' selection probability achievement over 10 runs using five-objectives engineering problem, Figure S2: Convergence progress of dominance solutions during evaluation process of engineering problem for 10 trials using Borg MOEA, Figure S3: Convergence progress of dominance solutions during evaluation process of engineering problem for 10 trials using ε-DSEA.