Developing Optimal Reservoir Operation for Multiple and Multipurpose Reservoirs Using Mathematical Programming

Over the last decades, the increasingwater demand has caused a number of problems, towhich reservoir operation optimization has been suggested as one of the best solutions. In this research, a model based onmixed integer linear programming (MILP) technique is developed for the systematic operation of multireservoirs that are used to cater for the different needs of the Tehran-Karaj plain. These reservoirs include Laar, Latian, and Karaj dams. The system configuration was accomplished through the nodes and arcs of the network flow model approach and system component implementation including sources, consumption, junctions, and the physical and hydraulic relationship between them. The following were performed via comprehensive developed software: system configuration, objective function and constraints formulation, linearization, determining penalty values, and setting priorities for each node and arc in the system. A comparison between the MILP developed model’s results against the periodic data shows 21.7% less overflow, 11.6% more outflow, and 15.9% more reservoir storage, respectively. The outcome of the MILP-based modeling indicates superior performance to the historical period.


Introduction
Water scarcity has become a serious concern in the process of urbanization and socioeconomic development, specifically in metropolis cities in which water resources are limited [1].Therefore, in order to control and optimize surface water resource utilization, several large dam and water transfer projects have been created or are being designed and implemented in most countries worldwide including Iran.Thus, it is crucial to pay particular attention to the projects' effective operation to obtain the utmost benefits and satisfaction from all the goals set earlier.
Uneven periodic and spatial distribution of fresh water and increasing population and social welfare have led to a rising need for water in Iran, principally in Tehran.
The limited resources of surface water supplied from the reservoirs of Laar, Latian, and Karaj are not sufficient to cover the water needs of this metropolitan city.Therefore in order to meet the demands, underground aquifers are significantly being used to supplement the water needs.However, during the arid years underground water usage increases drastically, which may have devastating effects on the quality and quantity of aquifers.Thus, it is important to ensure that the ground water is not overextracted, and this can probably be done by optimizing the use the surface water.A proper management of surface water may indirectly influence the use of scarce ground water.Therefore, the aim of this study is to develop a model that can lead to an optimal operation of the current reservoir for better water resource management.
Researchers have been continuously developing various methods of optimizing reservoir operation [2].Optimization processes depend on system features, data availability, objective function and constraint types, a number of constraints and variables, and physical relationships governing the system [3].Nonetheless, uneven periodic and spatial distribution of fresh water and increasing population and social welfare have led to a rising need for water in Iran, principally in Tehran.Therefore, one of the requirements in water resource management is to optimize the use of surface water resources.
Mixed integer linear programming (MILP) is an integer programming model that allows combining real, integer, and binary variables.It is considered an appropriate approach to formulate and solve problems [4].Since the majority of integer programming models have a similar structure to linear programming, all the available tools can be applied.Robustness, ease of expansion, and a multiobjective nature are additional key features of the mentioned method [5].In contrast, several other integer programming models are nonlinear (MINLP), and recently some researchers have applied this method for problem solving [6][7][8][9].
Since Laar, Latian, and Karaj reservoir systems are large scaled and contain many variables and nonlinear constrains during a long historical period, in order to optimize reservoir system in a specific area and to overcome the problems in linear programming and reservoir characteristics MILP is considered as an adequate approach.
In the present research, a MILP classical model was developed to operate the current reservoirs in the study area.The aim is to examine MILP's technical competence in deriving a public policy of reservoir operation.

Literature Review
Researchers have suggested optimization models to balance the conflicts among water users regarding allocation of water recourses [10].Various methods of optimization and their applications have been comprehensively assessed by numerous researchers such Simonovic [11] and Labadie [12]; issues addressed are water resource engineering and particularly reservoir operation.All techniques have advantages and disadvantages.Regarding the optimization of reservoir system modeling, among the optimal techniques is linear programming (LP) along with alternative methods.Significant advantages of LP method are its ability to solve large-scale problems in the best possible way, convergence towards the global optimum response, an initial solution being not necessary, easy to perform sensitivity analysis and uncomplicated problem solving [13].
Linear programming guarantees reaching optimum answers, but all equations (including objective function and constraints) must be linear [14].In reality, most functions are nonlinear in plenty of relevant cases in water resource management.Moreover, there is a problem to formulate linear programming for reservoir operation.
One of the drawbacks with LP formulation in reservoir operation and management is that reservoir storage continuity equation cannot explicitly control spillway, whereas perhaps when the reservoir is not full in optimized solutions some values are considered for the spillway [15].This issue also was reported by Shih and ReVelle [16].It is also possible that, in some reservoirs alongside the spillway outlet, the amount of seepage from the reservoir is considerable, and there is a nonlinear relationship between seepage and storage.
From all the progress made on LP we can point to binary linear programming, integer linear programming, and mixed integer linear programming, which are very useful to describe nonlinear and nonconvex statements in the objective function and the constraints [17].Many researchers like Needham et al. [18], Srinivasan et al. [19], Barros et al. [20], and Eslami and Qaderi [21] have used the MILP technique in water resource engineering.
Optimization models are usually used for hydroelectric power reservoir operations with various time scales [22].The mixed integer programming method has been widely implemented for planning and midterm scheduling [23], short-term scheduling [24], and real-time scheduling [25].In some research works, this technique is used to formulate and estimate water resource management planning under hydrological uncertainties [26,27].
Different integer programming means like the cut, branch, and bound methods can been applied to solve a range of problems in the field of water resource [28].But to practically solve difficult and large-scale optimization problems, software packages are used due to their high speed and accuracy [29].Ilog-Cplex model is one of the software packages with the capacity to solve integer programming as well as many linear or nonlinear programming problems.It also has the ability to interact with other models, such as GAMS, AMPL, and Microsoft Visual Studio [30].Researchers like Conejo et al. [31], García-González et al. [32], and Baslis and Bakirtzis [23] have utilized these packages to solve MILP and mixed integer nonlinear programming (MINLP).
Moraga et al. employed MILP approximation to solve a nonlinear model for planning midterm water transfer between two reservoirs [33].Tu et al. [34] developed a mixed integer nonlinear programming (MINLP) model to optimize hedging rules for reservoir operations.In other research works, El Mouatasim [35] used Boolean integer nonlinear programming (BINLP) for optimum pump performance in reservoir operation.Noory et al. [36] also successfully used the MILP model for optimum irrigation allocation and solving some multicrop planning issues.
The system related to reservoirs can be portrayed as a network of nodes and arcs.Nodes represent water storage points, diversion places, and intersections, while arcs represent the release of reservoirs, channel pipelines, and evaporation or other losses.This kind of optimization is based on the flow network model as used by various researchers in CalSim, Modsim, and OASIS software to model large systems [12,37,38].

Case Study
The study area covers the plain and city of Tehran, which are located on the southern slopes of the Alborz Mountains and include the basins and dam reservoirs of Karaj, Laar, and Latian.Other than Laar basin dam, the remaining study area is located on the southern slopes of the Alborz Mountains.Several dams are constructed on the rivers in Tehran to control the surface water flow and for supplying a part of urban drinking water and agriculture requirements.Three most important dams are Karaj, Laar, and Latian dams.Recently, Taleghan and Mamlu dams are built but due to the lack of historical data they have not been discussed here.Figure 1 shows the relationship between water resources (surface and ground water) and different uses in the Tehran plain.In this study, only surface water resources in operation are taken into account, including the Laar, Karaj, and Latian dam reservoirs.Table 1 lists the reservoirs' characteristics.Tehran's drinking water supplies are provided for by ground water and surface water sources.Ample resources are allocated yearly to meet drinking water, agriculture, industry, and green space needs.The extraction of ground water resources greatly increases throughout arid years, which has harmful effects on the quality and quantity of ground water resources.Since 1928, several plans have been implemented to supply the water needs of Tehran.Unfortunately, over time, with the increasing population and higher levels of various consumer needs, the implemented plans were unable to meet consumer needs.Several projects are currently under study and will be implemented to provide greater amounts of water to the city of Tehran from other areas.

Mathematical Formulation of the Problem
The mathematical formulation developed in this study is based on the MILP approach.It has been attempted to introduce a set of constraints that involves correct variables such that the spillway can be controlled and nonlinear equations are presented as a set of linear equations.Mixed integer linear programming (MILP) is applied to perform this action.

Mixed Integer Linear
Programming.The complexity of computing an integer programming problem depends on two factors: the number of integer variables and the problem structure.Presently, the most popular algorithm for solving integer programming problems is the branch and bound technique, whereby the key feature is the implied counting of current answers.Since the number of current solutions of an integer programming problem of integer number limit is countable, it seems natural for the counting method to be applied to obtain the optimal solution.Suppose that the relationship between the reservoir's release and storage is simulated as the curve fitting Figure 2.This curve must be linearized by the piecewise method.In Figure 2   is the initial storage,   is the beginning of the second part,   is the beginning of the third part, and   is the endpoint of the third part;   ,   , and   are the slopes of the lines in Sections 1, 2, and 3, respectively, and  1 ,   Indeed, the release-storage equation is a curve that is divided into several lines after the cut approximation, and each line has a separate equation: where   is the release of reservoir,   is the reservoir storage,  is the slope of the line, and ℎ 1 is the origin width.The release-storage equation comprises several lines.The overall relationship between release and storage is as follows: The relationship must be achieved so as to calculate the storage at any moment according to release.At any time in the release-storage relation, the storage of the reservoir must be in one section.Thus by using the integer variables, the above equation becomes where ℎ 1 , ℎ 2 , and ℎ 3 represent the intercept in parts 1, 2, and 3, and  1 ,  2 , and  3 are integer-independent variables, respectively, which can only have values of zero or one.The total reservoir storage is equal to the sum of the figures in each area of the reservoir: Also, reservoir storage should be located between the upper and lower bounds in every part: Integer variables can only have values of zero or one, and their sum must always be equal to one: According to the storage, one of the integer variables has a value of one and the other zero.So, the final constraints are defined as follows: where  is a constant value, depending on how a curve is approximated by the number of lines; the number of integer variables increases and the problem becomes more complicated.

System Configuration.
As previously mentioned, the optimization model of network flow is used to implement multiobjective optimization of problem reservoirs in the Tehran plain.System topology was executed based on this model, after which the whole system was solved based on MILP.
In this method, each system is shown as a combination of nodes and arcs.Nodes must represent sources, reservoirs, junctions, demands, and ground water, and the arcs represent the links between nodes potentially including channels (Figure 3).Some of the main justifications for using the network flow model for system configuration are the ease of expressing river basin systems as schematic figures, capability to use linear programming (and MILP), applicability to large-scale problems, and effortless ability to alter the system.In the network flow model, system components are implemented using nodes and arcs.Most of the analyses regarding reservoir systems are conducted via a type of network flow formulation which is called "Minimum Network Cost." Each node and arc have a minimum and maximum boundary and a value.In the network flow model, all nodes and channels enter into the objective function with a coefficient that reflects the priorities.
A sample water supply system is presented in Figure 3, which consists of nodes and arcs.In Figure 3 an overview of dams, refineries, tunnels, and channels of Tehran are demonstrated.Karaj dam is a multipurpose reservoir on the Karaj river.The intake of Bilaghan is located 23 km away from Karaj dam, which receives a part of water requirements from Karaj dam as in the river and Taleghan dam via pipes.Part of that is allocated to supply water of Tehran by pipeline and the other part is allocated to about 50 thousand hectares of agricultural land.
Laar dam is another important source of drinking water supply in Tehran.Laar dam also supplies water for agricultural irrigation in Mazandaran Province.Laar tunnel transfers outflow of Laar's dam to Latian reservoir which is approximately 140 million cubic meters.
Latian dam is located on Jajrood river of Tehran.After primary physical treatment, Latian release is transferred by two steel and concrete pipes to first and second refinery.Thus, it supplies about 30% of the drinking water of Tehran.Mamlu and Karaj channels are added in the second refinery and they supply drinking water of Tehran.

Constrains.
Constrains are defined as continuity at each node, continuity at reservoir, elevation-storage relations, and the amount of shortages in drinking water, industry, and agriculture.Constraints with the reservoir can be expressed as follows: The above equation represents the reservoir storage balance;  +1 is reservoir storage at the end of period ,   is reservoir storage at the beginning of period ,   is the total input in a period of , and   is the total release in a period of .
The second constraint regards the high boundary and low boundary of reservoir storage, which are defined to achieve suitable storage for flood control and dead storage, and so forth  ,min ≤   ≤  ,max for  = 1, . . ., , where   is reservoir storage in a period of ,  ,max is the maximum reservoir storage, and  ,min is the minimum reservoir storage in the same period.The third constraint is the supply of minimum downstream flow to maintain water quality, wildlife, and the environment: where   is release of the reservoir in a period of ,  ,max is the maximum release of the reservoir storage, and  ,min is the minimum release of the reservoir storage in the same period.It is possible in reservoirs that produce energy, for the relationships between the area-storage and elevationstorage are used to obtain the effective head.Water resource management depends on the water needs supply to a great extent.The equation of each necessary node for drinking water, industry, agriculture, and so forth is where  is a variable related to the input of a needed node and  is a variable related to the rate of shortage in the same node whose range must be defined according to the type of user needs;  is the amount of water needed.In addition to supplying needs, in each needed node the user must determine the penalty rate for shortages.The penalty amount depends on the type and quantity requirements.Penalties that are set for shortage of drinking water can be more than the penalties that are set for water shortage in agricultural or industrial purposes.According to the importance of each consumer, the penalty may vary.
The connections between nodes are placed by the arc.The arc represents the current between two nodes.Different types of arcs are designed in this model, which must be applied according to usage and type of connection.The considered constraint for each arc is as follows: where ℓ ,max is the maximum current flowing through the spillway, and ℓ is the actual amount of current passing through the spillway in a period .

Objective Function.
In this developed model objective function will be defined as Minimize ∑ ℓ  ℓ  ℓ in which "" is all of nodes and arcs that have a defined priority or penalty factor."" is amount of priority or penalty coefficient and "" stands for flow or shortage in each node and arc.The main objective in such issue is minimizing a series of designed phrases so penalties should be fitted to type of consumer and amount of consumption with a positive factor and besides amount of priorities should be defined with a negative factor.The greater the amount of penalties in each constraint, the less violation against constraints by the model.Also, the greater the value of priorities, the more attempts to satisfy the given constraints by the model.Therefore, user should define numeric amount of such priorities and penalties according to condition and also type and amount of consumption appropriately.The objective function is given as Minimize OBJ: where  is the number of reservoirs,  is a whole period of optimization, and   is the penalty coefficients or priority.
Clearly, the objective is to minimize the sum of a series of expressions, including controlled channels (CH), shortages of agricultural needs (HO), minimization of overflow (PL), maximization of reservoir storage (ST), and minimization of urban demand shortage (US). 1 is priority use of controlled channels,  2 is priority for agricultural purposes,  3 is penalty factor for overflow,  4 is priority to maintain reservoir storage, and  5 is penalty coefficient for the urban water supply shortages.The coefficients of these terms are determined by the designer.Penalties are greater if the deficiencies are more.Optimal release of each of the reservoirs should be allocated to various uses; of those, drinking water is the most important consumer of water in system of Tehran.The penalties and priorities have been set so that the need of Tehran drinking water is a priority.Therefore, the primary purpose of implementing the operation is minimization of water shortages in the study area as a way to increase the amount of the availability of water.The next priority of project is a high percentage of the available water supply in any time.So overflow rate is very important in assessing operating performance.Priorities are intended to minimize the amount of overflow and caused increased hydropower production and as a result, reservoir storage is maintained.
The MILP model presented in this paper is easily developable and very flexible, because penalty coefficients and objective priorities can be modified by the user.For example, we can consider more penalties for overflow than other objectives in order to minimize overflow as a primary objective.Inputs and outputs of the model are presented in the form of spreadsheet software that is easily transferable to other spreadsheets.
Presently a mixed integer linear programming model has been developed, which can be used for water resource system optimization in a desired range.This model is prepared based on the formulas given by ( 8) to (13) and the developed formulas' linearized relationship to the nonlinear equations ( 1) to (7) in this problem.The period is equal to 20 years to solve this problem; the number of integer variables is equal to 1440; the average number of variables is equal to 60960 and the number of constraints is 47040.

Results
In this study, we used the optimization model to simulate operation of multiple and multiobjective reservoir system, without specifying operation rules of reservoirs.It was done by determination of penalty weighting factors and prioritizing in the objective function.Therefore, we did not use the objective function to maximize economic efficiency; we used it to have a better simulation of reservoirs performance compared to historical data.In this study, the objective function is governing equation in simulation of multiple and multiobjective reservoirs operation.This equation consists of several decision variables including the amount of overflow, reservoirs storage, and the amount of water needed for drinking, agriculture, and industry in each period in each reservoir.For water resource system mode, the assumptions and criteria considered are drinking water shortage of 5% and between 10% and 30% shortage of agricultural water, and the return flow from drinking water consumption and agricultural water consumption are considered 60% and 25%, respectively.To maintain the reservoir water balance, water surface elevations at the beginning and end of the optimization period were consider as equal to historical period amount.The optimization is done monthly and a historical period of 20 years has been used for the model from 1983 to 2003.The calculated optimal release from the MILP model was compared against observations of the Laar, Latian, and Karaj reservoirs.The optimal release from Latian dam, Laar dam, and Karaj dam is provided in Figures 4, 5, and 6 respectively.
There is an almost good adaptation between the results of the MILP optimal release and the observation amount in dams of Karaj and Latian.The greatest differences in Karaj dam occur in the maximum peaks.Also in Latian dam, after 1995 the delay between the optimal and actual values can be seen, in which existence of periodic droughts can be the reason.A large change in Mazandaran agricultural needs is also one of the main reasons for the creation of the difference between the historical and optimal values in Laar dam.
To supply drinking water for Tehran, the average annual of 336 million cubic meters (MCM) (38.89%) is taken from Karaj dam, and 288 MCM (33.33%) is taken from Latian and Laar dams.The remaining needed water, that is, 240 MCM (27.78%), is provided by ground water resources (wells).The underground source is substantial to meet the demand of Tehran drinking water.These values are shown in Table 2.In comparison, Table 3 gives the calculated amount of water supplied by the Karaj, Latian and Laar, and the underground source.
The driest year in history was 2001, when totally an average of about 216 MCM of water was drawn from Karaj Dam and 224 MCM from Laar and Latian dam, while 390 MCM was sourced from the wells to supply the drinking water for Tehran as shown in Table 3.In comparison, according to the statistics provided by the Tehran Regional Water Company, while the amount withdrawn from the Karaj and Laar and Latian dam were 214 MCM and 224 MCM, respectively, the amounts withdrawn from the well were over 440 MCM.This amount differs significantly from the calculated one, which indicates over consumption of the ground water resources.
Figure 7 shows the comparison between the historical and calculated amount of water abstracted from the Karaj dam, Laar and Latian dam, and the underground source.It can be seen from this figure that the observed amount of water     abstracted from the wells is always higher than the calculated amount, whereas this is not the case for the Karaj dam and Laar and Latian dam. Figure 8 shows the amount of water abstracted from the underground source by historical and calculation (MILP) methods.Tables 4-6 show the comparison of the results obtained from the reservoirs' performance, reservoir storage, and overflow rates.
Table 4 shows that MILP model has more reservoir storage in all three reservoirs compared to historical period results.This increment of storage may lead to growth in hydropower generation.Table 5 demonstrates that the MILP model increases the outflow for each dam, whereby the maximum increment is 24.64% from the Karaj dam.The average long term of changes in overflows can be seen in Table 6; the results indicate that for Laar dam there is no overflow value in both historical period and MILP method.The negative sign in results means the reduction in overflow, in which in total around 21.7% reduction in all three dams is calculated.

Discussion and Conclusions
The historical data has shown that Kataj and Latian and Laar dams cannot supply the total water need of Tehran; thus the underground water sources are being used to answer all the demands.However, overusing the underground water has many negative effects quantitatively and qualitatively.In order to supply water needs and minimize the underground water consumption, a model based on MILP technique has been developed to achieve optimal operation of the mentioned reservoirs.In other words, by optimal operation of surface water we can achieve the minimum underground water utilization.The mentioned technique is developed based on the specifications of the existing reservoirs in the case study, data availability, and compatibility of the method with the operation problem type.The reservoir system's operating rules were prepared using weighting factors of penalties and priorities in the objective function.In this study, the maximum emphasis on weighting factors is on shortage of drinking water, the overflow of reservoirs, and

Figure 1 :
Figure 1: Schematic drawing of the relationship between resources and expenditures for Tehran plain.

Figure 2 :
Figure 2: The relationship between the release and storage of a dam's reservoir.

Figure 3 :
Figure 3: A schematic of the prepared water system (MILP model).

Figure 4 :
Figure 4: The optimal release and observations of Latian dam.

Figure 5 :
Figure 5: The optimal release and observations of Laar dam.

Figure 6 :
Figure 6: The optimal release and observations of Karaj dam.

Figure 7 :
Figure 7: The amount of water supplied from the dams of Karaj, Latian, Laar, and Tehran aquifer (in MCM).

Figure 8 :
Figure 8: Comparison between taken discharge from the aquifer by observations method and the calculation method (in MCM).

Table 1 :
Reservoir characteristics in the study area.

Table 2 :
The observed average annual amount of water abstracted from the dams of Karaj and Latian and Laar and wells (in MCM).

Table 3 :
The calculated average annual amount of water abstracted from the dams of Karaj and Latian and Laar and wells (in MCM).

Table 4 :
The comparison between the results of average storage of reservoir by using historical data and MILP method for Karaj, Latian, and Laar dams (in MCM).

Table 5 :
The comparison between the results of average outflow of reservoir by using historical data and MILP method for Karaj, Latian, and Laar dams (in cms).

Table 6 :
The comparison between the results of average overflow of reservoir by using historical data and MILP method for Karaj, Latian and Laar dams (in MCM).maintaining factors of reservoir storage level.The results obtained from the MILP model are more appropriate than historical operation results and it shows that the amount of water taken from aquifer is less in MILP calculation.The results of the developed MILP model weighed against the historical data show 21.7% less overflow, 11.6% more outflow, and 15.9% more reservoir storage.Having more reservoir storage leads to improvement in hydropower efficiency and increase in the stored water in the reservoirs. the