A Heuristic Dynamically Dimensioned Search with Sensitivity Information (hdds-s) and Application to River Basin Management

River basin simulation and multi-reservoir optimal operation have been critical for river basin management. Due to the intense interaction between human activities and river basin systems, the river basin model and multi-reservoir operation model are complicated with a large number of parameters. Therefore, fast and stable optimization algorithms are required for river basin management under the changing conditions of climate and current human activities. This study presents a new global optimization algorithm, named as heuristic dynamically dimensioned search with sensitivity information (HDDS-S), to effectively perform river basin simulation and multi-reservoir optimal operation during river basin management. The HDDS-S algorithm is built on the dynamically dimensioned search (DDS) algorithm; and has an improved computational efficiency while maintaining its search capacity compared to the original DDS algorithm. This is mainly due to the non-uniform probability assigned to each decision variable on the basis of its changing sensitivity to the optimization objectives during the adaptive change from global to local search with dimensionality reduced. This study evaluates the new algorithm by comparing its performance with the DDS algorithm on a river basin model calibration problem and a multi-reservoir optimal operation problem. The results obtained indicate that the HDDS-S algorithm outperforms the DDS algorithm in terms of search ability and computational efficiency in the two specific problems. In addition; similar to the DDS algorithm; the HDDS-S algorithm is easy to use as it does not require any parameter 2215 tuning and automatically adjusts its search to find good solutions given an available computational budget.


Introduction
It has been increasingly recognized that climate change and human activities have a profound impact on river basin management [1][2][3].The law of hydrologic cycle evolution, which could be explored by river basin simulation effectively, has a profound impact on river basin management in terms of efficiency of system optimization [4][5][6][7][8].Because the intense interactions between climate change and human activities lead to great changes in river quantity and quality regimes, river basin simulation becomes more difficult under the influence of many parameters.On the other hand, due to the uneven distribution of water resources in terms of time and space, a large number of reservoirs have been constructed in China as well as in many other countries.Although the joint operation of multi-reservoir systems could assist water resource optimal distribution in river basins, the complex topological structure of multiple reservoir systems in river basins makes the joint operation more difficult [9][10][11][12][13][14].
Because a large number of decision variables and constraints exist in the river basin model and the multi-reservoir operation model, developing high-performance optimization algorithms is critical for river basin model calibration problems and multi-reservoir optimal operation problems in river basin management.Use of conventional optimization algorithms to solve these optimization problems might result in slow convergence and low accuracy issues.Heuristic algorithms, which mimic human thinking, have become popular in recent years, as they are able to provide near-optimal solutions and can improve the speed of search [15,16].
There are many heuristic algorithms that are developed in recent years, such as Genetic Algorithms [17][18][19][20][21][22][23], Particle Swarm Optimization [24,25], Shuffled Complex Evolution (SCE) [26,27] and Ant Colony Optimization [28][29][30].The SCE algorithm is one of the popular optimization algorithms in the river basin model calibration over the past 10 years given that more than 300 different publications referenced the original SCE publications [26,27,31].However, various kinds of optimization algorithms need to utilize an extraordinary number of function (or model) evaluations to approximate the global optimum more accurately and it would considerable time to solve the complex optimization problems.
Generally, a wide variety of such problems could be solved in a simulation-optimization framework.The time required for one single simulation can vary from seconds to hours; the budget should be spent as efficiently as possible as it usually is limited.A global optimization algorithm called dynamically dimensioned search (DDS), which was proposed by Tolson and Shoemaker [32], is able to automatically refine and update optimal search space, and thus reduce the number of iterations required to reach satisfying solutions.Tolson and Shoemaker [32] compared the DDS and SCE algorithms from two aspects-The number of decision variables and the number of iterations, and found that the DDS algorithm is able to converge faster and get better solutions within a specified number of iterations, which proves its effectiveness for solving complex optimization problems.Due to the relatively good performance, simplicity and parsimony of the DDS algorithm, i.e., its ability to automatically scale the search space for good solutions within a specified computational time limit without requiring algorithm parameter adjustment, the DDS algorithm is ideally suited for solving computationally expensive optimization problems.
Most optimization algorithms treat each decision variable with an equal and constant probability to search and do not consider the variation of the probability in the search process, e.g., the SCE algorithm searches in the global space without dimension reduced.However, the DDS algorithm searches globally at the start of the search and conducts a more localized search as the number of iterations approaches the maximum allowable number of function evaluations, which is achieved by dynamically and probabilistically reducing the number of dimensions in the neighborhood and the reducing probability is distributed to each decision variable equally.Meanwhile, Hansen and Ostermeier [33] developed an adaptive evolutionary strategy to improve search efficiency, where decision variable co-variant matrix could be constructed and their automatic adjustment could be implemented in the search process.
This study presents a further developed DDS algorithm to solve high-dimensional optimization problems, i.e., river basin model calibration problem and multi-reservoir optimal operation problem in river basin management, with many parameters more efficiently.The newly developed global optimization algorithm, named as heuristic dynamically dimensioned search with sensitivity information (HDDS-S), is built on the DDS algorithm and could also search from global to local space with a reduced number of dimensions.However, during the adaptive change from global to local search with dimensionality reduced, instead of the reducing probability being distributed to each decision variable equally in the DDS algorithm, the non-uniform probability is assigned to each decision variable on the basis of its changing sensitivity to the optimization objectives in the HDDS-S algorithm.To evaluate the efficiency in solving high-dimensional and computationally expensive optimization problems in river basin management, the HDDS-S algorithm is tested extensively against the DDS algorithm on the river basin model calibration problem for the Tang-Wang River Basin in Jilin province, China, and the optimal operation problem for the DHF-GYG-SW (Dahuofang reservoir (DHF), Guanyinge reservoir (GYG), and Shenwo reservoir (SW)) multi-reservoir in the northeast of China.
This study presents the HDDS-S algorithm developed on the basis of the DDS algorithm to effectively perform river basin simulation and multi-reservoir optimal operation during river basin management.The original DDS algorithm is chosen as the benchmark optimization algorithm, that is, the HDDS-S algorithm is tested extensively against the DDS algorithm in solving high-dimensional and computationally expensive optimization problems in river basin management.The DDS algorithm was developed for the purpose of finding good global solutions, which are opposed to globally optimal solutions, within a specified number of maximum function evaluations.The key feature of the DDS algorithm was motivated by past experience from manual calibration of river basin models and reservoir operation models.A number of model parameters should be modified simultaneously early in the optimization exercise due to the fact that only relatively poor solutions could be obtained more often at this moment.However, as the optimization results improve, it becomes necessary to only modify a few parameters simultaneously so that the current gain in optimization results would not be lost.Additionally, through the observation of multi-dimensional decision variable optimization process, Tolson and Shoemaker [32] found that the number of decision variables that need to be optimized for improving the quality of solutions is decreasing as the values of decision variables approach their best values gradually.Therefore, unlike other heuristic algorithms that optimize each decision variable within a fixed space all the time, the DDS algorithm automatically sets the searching region size according to the current iteration count and user-specified maximum number of function (or model) evaluations to reflect the phenomenon that the number of decision variables to be optimized is decreasing, which is achieved by dynamically and probabilistically reducing the number of dimensions in the neighborhood and the reduced probability is distributed to each decision variable equally.The probability of choosing each decision variable to search based on a function of iteration count in the DDS algorithm could be calculated as: where ( , ) is the probability of choosing decision variable to search on iteration , is the total number of iterations.
On the other hand, in the original DDS algorithm, candidate solutions are created by perturbing the current solution values in the randomly selected dimensions only.These perturbation magnitudes are randomly sampled from a normal distribution with a mean of zero to keep the DDS algorithm as simple and parsimonious as possible.The DDS algorithm is a greedy type of algorithm since the current best solution is never updated with a solution that has an inferior value of the objective function.See Tolson and Shoemaker [32] for a detailed description of the DDS algorithm.

Heuristic Dynamically Dimensioned Search with Sensitivity Information
The regular used algorithms adopt equal probability search approach without dimension reduction, e.g., the SCE algorithm.The DDS algorithm could select dimensions to search randomly with a function of iteration count, and the number of selected dimensions is decreasing gradually with the increase of iteration count.Similar to the DDS algorithm, covariance matrix adaptation (CMA) and a multi-algorithm, genetically adaptive multi-objective (AMALGAM) algorithms could create the number of offspring points self-adaptively based on experiences, and all the points are searched with the same probabilities [33][34][35][36].
Heuristic dynamically dimensional search with sensitivity information (HDDS-S) improves the equal probability random search with dimension reduction in the DDS, CMA and AMALGAM algorithms, realizing the non-uniform probability search with the changing sensitivity of different dimensions to the objectives.The search approaches of different algorithms are shown in Figure 1.Because the number of decision variables that need to be optimized for improving the quality of solutions is decreasing as the values of decision variables approach their best values gradually, the changing process of decision variable sensitivity works on two principles: (1) the sensitivity of decision variable reduces with the increase of decision variables selected for searching on each iteration; (2) the sensitivity of decision variables is continuously cumulative with the increase of iteration number.
In the HDDS-S algorithm, to reflect the two principles, firstly, the sensitivity of each decision variable could be calculated as: where , is the sensitivity of decision variable ( = 1, ⋯ , ) on iteration ( = 1, ⋯ , ), is the number of decision variables selected for searching on iteration , is the collection of decision variables selected for searching on iteration , and are the total number of decision variables and iterations respectively.
Secondly, the cumulative sensitivity of decision variable could be calculated as: where , is the cumulative sensitivity of decision variable from 1 to iteration .
Similar to the DDS algorithm, the HDDS-S algorithm could search from global to local space with a reduced number of dimensions.However, instead of the reducing probability being distributed to each decision variable equally in the DDS algorithm, the non-uniform probability is assigned to each decision variable on the basis of its changing sensitivity to the optimization objectives in the HDDS-S algorithm.The probability of choosing a decision variable to search based on its cumulative sensitivity could be calculated as: where ( , + 1) is the probability of choosing decision variable to search on iteration + 1, and are the maximum and minimum cumulative sensitivity values among all decision variables from 0 to the iteration .The calculation of ( , + 1) reflects the dynamic adjustment of decision variable search strategy, that is, the non-uniform probability search with the changing sensitivity of different dimensions to the objectives.Therefore, the HDDS-S algorithm has the following characteristics: firstly, changing of the sensitivity matrix of multi-dimension parameters (decision variables) and cumulative sensitivity of each parameter are updated dynamically during the process of optimization according to the relationship between dimension selection and optimizing efficiency; secondly, the probability of choosing each parameter to optimize is calculated according to its cumulative sensitivity; finally, the optimization strategy of choosing parameters with non-uniform probability during the process of parameter optimization is realized to improve convergence rate and optimization efficiency.Additionally, similar to the DDS algorithm, the only algorithm parameter to set in the HDDS-S algorithm is the scalar neighbourhood size perturbation that defines the random perturbation size standard deviation as a fraction of the decision variable range.Tolson and Shoemaker [32] recommended the value of r parameter as 0.2 (and used in this study) because this yields a sampling range that practically spans the normalized decision variable range for a current decision variable value halfway between the minimum and maximum.
The complete HDDS-S algorithm pseudo-code for minimization can be seen in the Appendix Pseudo-code of HDDS-S.

Model Description and Study Area
Almost all river basin models contain effective physical and/or conceptual model parameters that are either difficult or impossible to directly measure.Optimization algorithms could be used to automate the calibration, i.e., automatic calibration is defined as an optimization algorithm-based search for a set of river basin model parameter values that minimize the model prediction errors relative to available measured data for the system being modeled process.Tolson and Shoemaker [32] compared the performance of the DDS algorithm to that of the SCE algorithm, which is currently one of the most commonly applied algorithms for automatic calibration of river basin models, with multiple optimization test functions as well as real and synthetic SWAT2000 model automatic calibration formulations.In that study, numerical results demonstrate that the DDS algorithm is a more computationally efficient and robust optimization algorithm than the SCE algorithm in the context of river basin model automatic calibration.Therefore, the HDDS-S algorithm developed on the basis of DDS is evaluated by comparing its performance with the DDS algorithm on a river basin model calibration problem in this paper.
Meanwhile, river basin simulation is closely related to multi-reservoir optimal operation.Firstly, the future runoff series under the changing conditions could be generated on the basis of the law of hydrologic cycle evolution, which could be explored with river basin model effectively, e.g., SWAT model.Secondly, the future runoff series need to be used as the inputs of the multi-reservoir operation model to explore the impacts of changing conditions on multi-reservoir operation, determine the operation rules, and assist in making multi-reservoir optimal operation decisions under changing conditions.As with river basin simulation, multi-reservoir optimal operation has been critical for river basin management; the multi-reservoir operation model contains large number of decision variables to be calibrated, i.e., large number of the curves and the number of time periods, which increases the complexity and problem, and poses a significant challenge for most optimization tools currently available [9][10][11][12][13][14]. To evaluate the universality of the HDDS-S algorithm in river basin management, the HDDS-S algorithm is also evaluated by comparing its performance with the DDS algorithm on a multi-reservoir optimal operation problem in this paper.
Two cases are used to extensively test our new HDDS-S algorithm against the DDS algorithm in this paper: (1) the river basin model calibration for Tang-Wang River Basin; and (2) the optimal operation for the DHF-GYG-SW multi-reservoir.The main difference between the multi-reservoir operation model and river basin model is that the decision variables in multi-reservoir operation models are reservoir operation rule curves, which consist of a series of storage volumes or levels at different periods, and the aim of multi-reservoir optimal operation is to establish reservoir operation rule curves at the planning/design stage and provide long-term operation guidelines for reservoir managers to meet expected water demands [37][38][39].Therefore, the decision variables, which will be optimized in the two optimization problems, are parameters of the river basin model and multi-reservoir optimal operation model, and the dimension being varied is the number of model parameter values being changed to generate a new search neighbourhood.

Overview of the SWAT Model
The SWAT model is a catchment-scale semi-distributed river basin model developed by the Agricultural Research Service of the United States Department of Agriculture; its detailed description can be seen in many references, such as [40,41].The SWAT model can be run at a daily time step when the Soil Conservation Service (SCS) curve number method is used to calculate surface runoff, and it is not designed to accurately simulate detailed, single-event flood routing.In this paper, the river basin model calibration problem is formulated based on the set of the SWAT2005 flow-related parameters, which impact snowmelt, surface runoff, groundwater, and lateral flow simulation.
Because SWAT is a semi-distributed river basin model and Hydrological Responding Units (HRUs) divided with geographic information are its basic units, some parameters are not regarded as spatial variables but rather constant values across all spatial units and some other parameters are spatially varied among different HRUs.For example, those parameters related to snowmelt have constant values in all spatial units, and could be calibrated one by one.However, some parameters, such as SCS curve numbers, could be assigned different values for different spatial units, which significantly increase the complexity and computational model requirement of parameter calibration.In this paper, a single factor is used to represent spatial variation, by increasing or decreasing spatially varying parameter values from their base or default values.For spatially varied parameter, this approach maintains the relative differences in the base or default parameter values assigned to different spatial units.
The flow-related parameters of the SWAT model and their ranges are listed in Table 1, where the parameter ranges are based mainly on the default ranges in the SWAT2005 model documentation.

Evaluation Criterion and Objective Function
The mean relative error ( ), the coefficient of determination ( ), and the Nash-Sutcliffe efficiency ( ) are used to evaluate the simulated streamflows with the observed streamflows.is computed according to Equation (6): where and are the means of the simulated streamflows and the observed streamflows, respectively." values of 0 indicate a perfect fit.Positive values indicate model overestimation bias, and negative values indicate model underestimation bias [42]." is computed according to Equation ( 7): where is the th simulated value for the streamflows; is the th observation for the streamflows; and is the total number of observations." describes the portion of the variance in measured data explained by the model.ranges from 0 to 1, with lower values indicating more error variance, and typically = 1 is considered the optimal value [43]." is computed according to Equation (8): The ranges between −∞ and 1.0, with = 1.0 as the optimal value.According to Luo et al. [44], values between 0.0 and 1.0 are generally viewed as acceptable levels of performance, whereas a value less than 0.0 means that the mean observed value is a better predictor than the simulated value, which indicates unacceptable performance.
To evaluate the HDDS-S algorithm by comparing its performance with that of the DDS algorithm on river basin model calibration, the sum of squared differences between the simulated streamflows and the observed streamflows is used as the objective function, which is shown in Equation ( 9): where is the sum of squared differences between the simulated streamflows and the observed stream flows.

Multi-Reservoir Optimal Operation Model
Multi-Reservoir Operation Rule Multi-reservoir operation rule consists of operation rule curves and rationing factors for each water demand.Details of the operation rule curves and their corresponding water supply rationing factors are illustrated in Figure 2 and Table 2.The forms of reservoir operation rule curves based on 36 10-day periods are shown schematically in Figure 2. The operation rule curves are defined according to reservoir storage, that is, the dynamic water storage of the reservoir is taken as the single, influential factor for reservoir operation.In this paper, the active water storage of reservoir is divided into three parts, zone I, zone II and zone III, by two operation rule curves for industrial (D1) and agricultural (D2) demands.Each of the water demands has a related rule curve and a rationing factor.Because different kinds of water demands require different reliabilities and different degrees of priority in practice, the rationing factors usually have different values.For example, assume water rationing factor is α1 for the operation rule curve for industrial demand (D1) in Figure 2. When the water storage of the reservoir lies in zone I or zone II, water demand D1 is fully met, i.e., the amount of water supply is D1; when the water storage of reservoir is in zone III, the amount of water supply is calculated by applying the water rationing factor, i.e., α1 × D1.

Multi-Reservoir Operation Model
The multi-reservoir optimal operation problem is formulated as a single-objective optimization problem, where decision variables are the water storages at different periods on different operation curves, including industrial and agricultural curves of XN-2 (the virtual reservoir which GYG and SW reservoirs are aggregated into) and DHF reservoirs.Meanwhile, the operation for the entire system should meet the guaranteed rate of each water demand, i.e., the guaranteed water-supply rate of industry and agriculture are above 95% and 50%, respectively.The guaranteed rate of each water demand can be calculated by statistics from the operation result, shown in Equation ( 10): where is the vector of decision variables, i.e., the water storages at different periods on an operation rule curve; is the guaranteed rate for water demand (industrial water demand when = 1, agricultural water demand when = 2), which measures the frequency of water demand being satisfied during entire time periods and is used as an indicator to reflect water supply capacity of the system; is the entire operation period for demand (industrial demand when = 1, agricultural demand when = 2); similar to , is the time period when water demand is being satisfied.Additionally, the number of depth damages, which indicate that the reservoir water storage has been lowered to dead water storage, should be as few as possible during the operation for the entire system.Similarly, the number of depth damages could be calculated by statistics from the operation result.
In this paper, the objectives are set to minimize the number of depth damages, ans improve the guaranteed rate of each water demand during the operation for the entire water supply system: where , is the guaranteed rate for water demand in subsystem (Sub A when = 1, Sub B when = 2, and Sub C when = 3); similarly, , is the number of depth damages for water demand in subsystem .
For the multi-reservoir operation optimization problem, the mass balance equations are: min ≤ ≤ max , min ≤ ≤ max (14) where is the initial water storage at the beginning of period ; is the ending water storage at the end of period ; , , , and are inflow, delivery for water use, spill and evapotranspiration loss, respectively; and max and min are the maximum and minimum storage, respectively.

Tang-Wang River Basin
The Tang-Wang River is the first level tributary of the left bank of Song-Hua River, China.The total length of the Tang-Wang River is approximately 509 km.Its basin drains an area of 20,383 km 2 .The study area is limited to the upper-middle stream region of the Tang-Wang River above Yixin hydrologic station, herein still referred to as Tang-Wang River Basin.The river basin drains an area of 11,376 km 2 .The basin boundary and associated SWAT model sub-basin boundaries are presented in Figure 3.There are 87 sub-basins defined in the Tang-Wang River Basin, where 16 rain gauges and 1 streamflow gauge are located.Daily Meteorologic data (temperature, solar radiation, weed speed, relative humidity), daily precipitation data and stream flow data used in SWAT Model Calibration for Tang-Wang River Basin are listed in Table 3.To evaluate the HDDS-S algorithm by comparing its performance with that of the original DDS algorithm on river basin model calibration, two scenarios are constructed: (1) model auto-calibration with the HDDS-S algorithm; and the (2) model auto-calibration with the DDS algorithm.The Yixin streamflow gauge was chosen as calibration station, and streamflow data from the 1979-2001 period were used to calibrate the model parameters.For each scenario, the first two months (January and February) are used as a warm up period for model simulation.The rest of the time periods in each scenario are used to assess the model's performance in the parameter calibration process.

DHF-GYG-SW Multi-Reservoir
This paper use thes DHF-GYG-SW multi-reservoir located in Huntai Basin as a case study, which consists of the Dahuofang reservoir (DHF), Guanyinge reservoir (GYG), and Shenwo reservoir (SW).Huntai basin, including Hun River and Taizi River, is located in a rapidly urbanizing region of Northeast China (Figure 4).DHF is a carry-over storage and is located in the trunk Stream of Hun River.It drains an area of 5437 km 2 , and the total length of the Hun River is approximately 169 km within the basin.GYG and SW are cascade reservoirs and located in the trunk stream of the Taizi River.GYG is a carry-over storage and drains an area of 2795 km 2 , while SW drains an area of 6175 km 2 with incomplete year regulation performance.The reservoir characteristics, annual average inflow and water supply tasks are illustrated in Table 4.The whole multi-reservoir system consists of three parts: subsystem A (Sub-A), which contains DHF and DHF-Sanchahe (SCH) interval; subsystem B (Sub-B), which contains GYG, GYG-SW interval, SW reservoir, and SW-SCH interval; and subsystem C (Sub-C), which contains SCH-Daliaohekou (DLHK) interval.The corresponding topological structure is generated, as shown in Figure 5a.Sub-A is supplied with water from DHF and uses DHF operation rules to decide how to supply water to each user in this subsystem.Sub-B is supplied with water from GYG and SW, which constitute cascade reservoirs and so are aggregated into one virtual reservoir (XN-2 reservoir), as shown in Figure 5b.Sub-C is supplied with water from the parallel reservoir, i.e., XN-2 and DHF, the common tasks in Sub-C should be allocated to XN-2 and DHF.The water ratio method is employed to allocate the common tasks in Sub-C, which is a simple and efficient method shown in Figure 5c,d.The calculation of water division coefficient α is shown in Equations ( 15)-( 17): where , is the initial water storage of DHF at the beginning of period ; , and , are the inflow into DHF and water supply from DHF during the period , respectively; , is the initial water storage of virtual aggregate reservoir XN-2 at the beginning of period , i.e., the sum water storage of GYG and SW; , and , are the inflow into XN-2 and water supply from XN-2 during the period , i.e., the sum inflow into GYG and SW, the sum water supply from GYG and SW, respectively.Under the same number of function evaluations ( = 1000), the average cumulative varied times for each parameter during model auto-calibration process in different scenarios through 30 optimization trials are shown in Figure 9, where the average cumulative varied times in S0 indicate the cumulative sensitivity of each parameter.= 1000) in different scenarios through 30 optimization trials, where the description of each parameter in the x axis could be found in Table 1.

River Basin Model Calibration
The average cumulative varied times for each parameter during model auto-calibration process under the same number of function evaluations ( = 1000) through 30 optimization trials are different between the two scenarios, i.e., the DDS algorithm, which selects dimension for perturbation completely at random with a uniform probability, could not reflect the sensitivity information of each parameter.The details are indicated from two aspects.
(1) As indicated in the cumulative varied times in the case of the HDDS-S algorithm, some parameters need to be varied many times due to strong sensitivities, e.g., SFTMP, ALPHA_BF and ESCO.However, the cumulative varied times of these parameters in the DDS algorithm are less.
(2) Similarly, some parameters need to be varied fewer times due to weak sensitivities, e.g., SMTMP, SMFMX, SURLAG and CN2.However, their cumulative varied times in the DDS algorithm are greater than.
Under the same number of function evaluations ( = 1000) through 30 optimization trials, average calibration processes of model parameters SFTMP, SMFMX, SURLAG and ESCO, in different scenarios are shown in Figure 10.It is indicated from Figures 9 and 10 that the value of the most sensitive parameter SFTMP needs to be varied through the whole calibration process, while its calibration is stopped when reaching some iteration number in the DDS algorithm.Similarly, the value of some more sensitive parameters, such as ESCO also needs to be varied many times, while it is calibrated with fewer variation times, scatter variation moments and obvious stagnation phenomenon in the DDS algorithm.However, less sensitive parameters SURLAG and SMFMX needs to be varied fewer times, while it is calibrated with more variation times, scatter variation moments and obvious stagnation phenomenon in the DDS algorithm.
Therefore, no matter how sensitive these parameters are, they could be calibrated to their optimal values much faster with the HDDS-S algorithm than the DDS algorithm.
River basin model calibration results show that the computation time and efficiency could be improved in the HDDS-S algorithm, including the enhanced convergence rate of objective function, concentrated varied moments, reduced blind variation and stagnation phenomenon, compared to the DDS algorithm.Therefore, the HDDS-S algorithm is more efficient than the DDS algorithm in the river basin model calibration for the specific catchment.

Multi-Reservoir Optimal Operation
The optimal operation problem for DHF-GYG-SW multi-reservoir is used to evaluate the universality of the HDDS-S algorithm in river basin management.The input data are inflow and water demand data, including the industrial and agricultural water demands in the future planning year, i.e., 2030, and the history inflow from 1956 to 2006.Each calculation year is divided into 24 time periods (with ten days as scheduling horizon from April to September, and a month as scheduling horizon in the remaining months).Since the industrial water supply occurs through the whole year, there are twenty-four decision variables for industrial water supply.The agricultural water supply occurs only during the periods from the second ten-day of April to the first ten-day of September, thus there are fifteen decision variables for agricultural water supply.Therefore, the entire operation system has four rule curves and (24 + 15) × 2 = 78 decision variables in total.
Similarly, based on the two scenarios S0 and S1, i.e., multi-reservoir optimal operation in the case of the HDDS-S algorithm and multi-reservoir optimal operation in the DDS algorithm respectively, average performances of the two algorithms in multi-reservoir optimal operation are evaluated under the same iteration number ( = 10,000) through 30 optimization trials from three aspects, including average required computation time, average objective function value, and average guaranteed water-supply rate on Intel Core™ 2 Duo, with a 2.66-GHz CPU.
First of all, the rationality of operation rules derived after 10,000 function evaluations in one random optimization trial in the case of the HDDS-S algorithm is evaluated, and the joint scheduling charts are presented in Figure 11.
As industrial water demand is more urgent than agricultural water demand, and the agricultural water would be limited primarily when reservoir storage is insufficient, the water supply rule curve for agriculture is above the operation rule curve for industry, as shown in Figure 10.In the scheduling cycle duration, each rule curve decreases before the flood season and then has a different degree of uplift after the flood season.In order to prepare for the large inflow of the summer flood season and to reduce surplus water, the reservoir active storage before the flood season should be lowered and then the number of limited water supply is reduced, i.e., the operation rule curves are decreased before the flood season.After the flood season, when the inflow decreases gradually, in order to ensure that no depth damage happens in dry years, the curves are raised to fill reservoirs as early as possible and the number of limited water supply is increased.In addition, the largest agricultural water use happens from April to May, thus the operation rule curve of agriculture is higher in the duration to conserve water and mitigate potential future shortages for industrial demand in some dry years.Specifically, referred to DHF reservoir, its water supply rule curve for industry falls in the duration from March to April due to spring floods produced by melting snow, meanwhile, its water supply rule curve for industry raises in May because the largest agricultural water use happens from April to May and much water has been supplied to agriculture.(2) In Sub B, the guaranteed rate for industrial and agricultural water demands are 95.81% and 75.46% respectively; (3) In Sub C, the guaranteed rate for industrial and agricultural water demands are 96.08% and 62.56% respectively.The average guaranteed rates in all three subsystems meet the requirement, i.e., guaranteed rate for industrial and agricultural water demands exceed 95% and 50% respectively.
Therefore, as mentioned above, the derived operational rule curves of DHF reservoir and XN-2 reservoir in the HDDS-S algorithm and the scheduling results are reasonable.
In different scenarios, under the same number of function evaluations ( = 10,000), average required computation time and objective function value through 30 random optimization trials are shown in Table 5.Under the same number of function evaluations ( = 10,000), though larger differences of average required computation time exist in different scenarios, average objective function value decreases from 25.68 in the DDS algorithm to 9.13 in the HDDS-S algorithm, which improve obviously.
Average guaranteed water-supply rates for all water supply area under the same number of function evaluations ( = 10,000) through 30 random optimization trials in different scenarios are shown in Table 6.Under the same number of function evaluations ( = 10,000), though few differences of average guaranteed rate for industrial water demand exist in different scenarios, i.e., the guaranteed rate for industrial water demand in all three subsystems exceed 95%, the average guaranteed rates for agricultural water demand in all three subsystems increase from 67.52%, 64.78% and 50.94% in the DDS algorithm to 75.35%, 75.46% and 62.54% in the HDDS-S algorithm in Sub A, Sub B and Sub C respectively, which improved obviously.Optimal operation of DHF-GYG-SW multi-reservoir results show that the derived operational rule curves of DHF reservoir and XN-2 reservoir and the scheduling results from one random optimization trial are reasonable in the case of the HDDS-S algorithm, while average guaranteed water-supply rates in all subsystems could be improved in the HDDS-S algorithm under the same number of function evaluations ( = 10,000) through 30 random optimization trials compared to the DDS algorithm.Therefore, the HDDS-S algorithm is able to derive good solutions and thus is effective in the specific large multi-reservoir optimal operation.

Discussions and Conclusions
The DDS algorithm could automatically refine and update optimal search space, i.e., searches globally at the start of the search and conducts a more localized search as the number of iterations approaches the maximum allowable number of function evaluations, and automatically adjust and exhibit good performance within the user's timeframe for calibration to obtain a quick and approximate result, which is achieved by dynamically and probabilistically reducing the number of dimensions in the neighborhood and the reducing probability is distributed to each decision variable equally [45,46].With multiple optimization test functions as well as real and synthetic SWAT2000 model automatic calibration formulations, Tolson and Shoemaker [32] have proved that the DDS algorithm could return good solutions faster than the SCE algorithm within the maximum allowable number of function evaluations.
In this paper, the DDS algorithm is improved with reference to the changing sensitivity information of multi-dimension decision variables, and the HDDS-S algorithm is proposed based on the sensitivity analysis.The HDDS-S algorithm consists of three steps: firstly, the changing sensitivity matrix of multi-dimension parameters and cumulative sensitivity of each parameter are updated dynamically during the process of optimization according to the relationship between dimension selection and optimizing efficiency; secondly, the probability of choosing each parameter is calculated during the process of optimization according to its cumulative sensitivity; finally, the optimization strategy with non-uniform probability of choosing parameters during the process of parameter optimization is realized to improve convergence rate and optimizing efficiency.The HDDS-S algorithm is tested extensively against the DDS algorithm on the river basin model calibration problem for Tang-Wang River Basin in Jilin province, China, and the optimal operation problem for DHF-GYG-SW multi-reservoir in the northeast of China.
In the river basin model calibration problem for Tang-Wang River Basin, although there is little difference of model performances between the DDS algorithm and the HDDS-S algorithm because the simulation results will approximate their optimal values after 1000 function evaluations to a large extent, the average required computation time reduces from 115.74 h ( = 796) in the case of the DDS algorithm, to 69.90 h ( = 408) in the HDDS-S algorithm, through 30 optimization trials under the same simulation accuracy ( = 370,000).Meanwhile, due to the consideration of the changing sensitivity information of each dimension decision variable, under the same iteration number, the convergence rate of the objective function is improved in the HDDS-S algorithm compared to the DDS algorithm, and there are few differences of the HDDS-S algorithm convergence processes among different optimization trials, which indicate the reliability and stability of the HDDS-S algorithm in solving the model auto-calibration problem.Additionally, in the optimal operation problem for the DHF-GYG-SW multi-reservoir, under the same number of function evaluations ( = 10,000), though few differences of average guaranteed rate for industrial water demand exist in different scenarios, i.e., the guaranteed rate for industrial water demand in all three subsystems exceed 95%, the average guaranteed rates for agricultural water demand in all three subsystems increase from 67.52%, 64.78% and 50.94% in the DDS algorithm to 75.35%, 75.46% and 62.54% in the HDDS-S algorithm in Sub A, Sub B and Sub C respectively.Overall, results show that the HDDS-S algorithm could be more efficient and effective than the DDS algorithm in the two specific problems.
Since one algorithm may be better for calibration of the river basin model for only one catchment or one specific multi-reservoir optimal operation, we will choose other catchments and multi-reservoir operation systems with different characteristics including climate, topography and human activities, to examine the performance of the HDDS-S algorithm under different conditions in future studies.and Yu Li provided their support and guidance throughout the manuscript writing.Their timely suggestions have greatly improved the manuscript.All authors reviewed the final manuscript.

Pseudo-Code of HDDS-S
The complete HDDS-S algorithm pseudo-code for minimization is provided in the following: STEP 3. Referred to all decision variables, according to their cumulative sensitivities from beginning, select of the decision variables to search, i.e., to be included in neighborhood, , take decision variable ( = 1, ⋯ , ) as an example to illustrate the process:

Figure 1 .
Figure 1.The search approaches of different algorithms.(a) Equal probability search approach without dimension reduction; (b) Equal probability search approach with dimension reduction; (c) Non-uniform probability search with dimension reduction.

Figure 4 .
Figure 4.The location of the adopted water supply multi-reservoir.

Figure 5 .
Figure 5. Reservoir aggregation and common water supply task allocation, (a) initial topological structure; (b) formation of virtual aggregate reservoir; (c) final topological structure after reservoir aggregation; (d) common water supply task allocation in the final topological structure.
The river basin model calibration problem for the Tang-Wang River Basin is used to validate effectiveness of the new HDDS-S algorithm.The Yixin hydrologic station, located at the outlet of Tang-Wang River Basin, is chosen as the calibration station, and the streamflow data from the 1979-2001 periods are used to calibrate model parameters.Based on the two scenarios S0 and S1, i.e., model auto-calibration with the HDDS-S and DDS algorithms respectively, performances of the two algorithms in the river basin model parameter auto-calibration are evaluated on Intel Core™ 2 Duo, with a 2.66-GHz CPU.In different scenarios, average required computation time under the same number of function evaluations, i.e., = 1000 , or the same simulation accuracy, i.e., = 370,000 , through 30 optimization trials are shown in Figure6.

Figure 6 .
Figure 6.Average required computation time under the same number of function evaluations ( = 1000) and the same simulation accuracy ( = 370,000) in different scenarios through 30 optimization trials.

Figure 7 .
Figure 7. Observed and average simulated monthly streamflows during the period from 1979 to 2001 at the Yixin hydrologic station through 30 optimization trials.Seen from Figures 6 and 7, under the same number of function evaluations ( = 1000), the average required computation time increases from 145.40 h in the case of the DDS algorithm to 171.32 h in the HDDS-S algorithm.Although the model performance improves from 11.79% ( ) in the case of the DDS algorithm to 9.94% in the HDDS-S algorithm at the Yixin hydrologic station, there is little difference of model performances between the DDS algorithm and the HDDS-S algorithm in terms of the other two evaluation criterions because the simulation results will approximate their optimal values after 1000 function evaluations to a large extent.Under the same simulation accuracy ( = 370,000), the average required computation time reduces from 115.74 h ( = 796) in the case of the DDS algorithm to 69.90 h ( = 408) in the HDDS-S algorithm through 30 optimization trials.Although the average computation time in the case of the HDDS-S algorithm increases due to the additive sensitivity computation of each dimension decision variable compared to the DDS algorithm under the same iteration number, the average computation time in the case of the HDDS-S algorithm reduces greatly compared to the DDS algorithm under the same simulation accuracy.The performances of the two algorithms in river basin model calibration are evaluated not only from average required computation time under the same number of function evaluations, average required computation time under the same simulation accuracy, and average simulation results at Yixin hydrologic station under the same number of function evaluations, but also from average convergence process of the objective function and average calibration process of model parameters under the same number of function evaluations.Convergence processes of the objective function under the same number of function evaluations ( = 1000) through 30 optimization trials in different scenarios are shown in Figure 8.It is further indicated that due to the consideration of the changing sensitivity information of each dimension decision variable, under the same number of function evaluations, convergence rate of

Figure 8 .
Figure 8. Convergence processes of the objective function under the same number of function evaluations ( = 1000) in different scenarios through 30 optimization trials, (a) and (b) convergence processes in two random optimization trials; (c) average convergence process through 30 optimization trials.

Figure 9 .
Figure 9. Average cumulative varied times for each parameter during model auto-calibration process under the same number of function evaluations (= 1000) in different scenarios through 30 optimization trials, where the description of each parameter in the x axis could be found in Table1.

Figure 11 .
Figure 11.The derived operational rule curves with the HDDS-S algorithm after 10,000 function evaluations in one random optimization trial, (a) DHF reservoir operation rule curve; (b) XN-2 reservoir operation rule curve.
a Parameters are multiplicative factors used to adjust the spatial variation across all model units.

Table 2 .
Operation rule implied by rule curves.

Table 3 .
Hydrometeorological data for Tang-Wang River Basin.

Table 4 .
Reservoir characteristics and tasks of water supply.

Table 5 .
Average computation time and objective function value under the same iteration number ( = 10,000) through 30 random optimization trials in different scenarios.

Table 6 .
Average guaranteed water-supply rate under the same iteration number ( = 10,000) through 30 random optimization trials in different scenarios.
Calculate the sensitivity on iteration : , . Calculate the cumulative sensitivity from beginning to iteration : , . Calculate the probability of choosing to search based on its cumulative sensitivity Update objective function iteration counter, = + 1, and check stopping criterion:  If = , STOP, save (e.g., & ). Else, go to STEP 3.