Advances in using mathematical optimization to manage floods with assessment of possible benefits using a case study

The use of computer models is a common component in the process of creating river basin management plans. Early review papers of the available models and their capabilities have been provided by Yeh (1985) and Labadie (2004), while more recent reviews are available from Rani & Moreira (2010) and Dobson et al. (2019), which focus on the advantages of optimization models over those based on simulation. In spite of the long history of diverse optimization model development, their actual use in real-time reservoir operation is almost non-existent. Some of the principal reasons for this are (a) the lack of reliable runoff forecast and (b) the inability of most optimization models to combine reservoir release decisions with their routed effects at downstream control points. The effects of channel routing should be combined with the effects of additional tributaries and downstream hydraulic structures. Also, the propagation of reservoir releases to the downstream control points may take longer than the short length of the calculation time steps that are required for adequate modeling of floods, thus requiring simultaneous solutions over multiple calculation time steps – a feature that is not available to simulation models that have traditionally been used to model flood events. This paper deals with the issue (b), i.e., the focus of the paper is not the runoff forecasting capabilities, but rather the required technical capabilities of river basin management models, where reservoir operation takes the center stage as the principal mechanism for taming the otherwise often hostile natural flow regime. Since reservoir operation is an essential part of river basin management, the terms ‘reservoir operation models’ and ‘river basin management models’ in this paper are assumed to refer to the same class of mathematical tools.

The first reservoir operation models such as the HEC-5 model (US Corps of Engineers 1998) appeared in the literature more than 50 years ago, and most of them were initially based on the concept of using simulation tools that allowed for inclusion of ‘what-if’ rules. Although some simulation models such as the Mike Hydro Basin (DHI 2023) and RIBASIM (DeltaRes 2023) are still used, their use has become increasingly difficult with the growth of the complexity of river basin management issues. These models often fail to simultaneously apply the zoning management concept on multiple reservoirs and multiple water users, so as to include not only the storage management policies, but also water rationing policies and deficit sharing in dry years as a function of the available storage and hydrologic conditions in the basin. Efforts to offset these difficulties led to the development of models based on the Network Flow Algorithms (NFAs), which are also based on the theory of linear programming (LP). The first NFA-based model was launched by the Texas Water Development Board (SIMYIELD) in the 1970s, followed by a number of models that have become popular over time, such as the MODSIM (Labadie et al. 1986), Aquatool (Andreu et al. 1996), E-Water (EWater Academy 2022), or the REALM model (Victoria State Government 2022). These models all have the capability to apply the zoning concept with reservoir rule curves in combination with the zoning concept on demands, such that demand management is also coupled with the management of water supply. Typically, the purpose of a reservoir operation model is to allocate water based on the allocation policies that are provided by the model users. While this approach may be helpful for managing storage during droughts, its usefulness is questionable during floods, where several upstream reservoirs can be managed concurrently, while inflow hydrographs into each reservoir always have unique temporal and spatial distributions. An additional limitation of NFAs is that most of them do not allow losses or gains of flow along a flow link in the network, which significantly complicates the inclusion of hydrological channel routing in the solution process.

Some river basin models, such as the WEAP model (Yates et al. 2005), also include an internal rainfall-runoff module that can generate runoff estimates. However, when it comes to river basin management models, this feature is more of an exception rather than the rule. Most river basin management models do not generate runoff estimates based on rainfall data, but rather rely on a previously prepared runoff time series for a period of multiple years or for individual hydrologic events. This makes the evaluation of the reservoir operation models easier, since the model performance can be separated from the ability to accurately predict runoff forecasts. The rainfall-runoff models should be calibrated and validated independently, and they are not discussed further in this paper. The principal question that this paper tries to address is the ability to improve real-time flood management assuming that an accurate short-term (3 days) runoff forecast is available, and to demonstrate the fact that an optimization-based computer model can be a powerful tool for real-time decision-making if coupled with a reliable runoff forecasting tool.

The use of optimization algorithms for reservoir operation during floods was first postulated as a conceptual LP model by Windsor (1973), and later refined by Needham et al. (2000). These applications were based on the development of flood damage functions, which were piece-wise linearized and provided to the LP optimization routines that are part of the model. Both papers detail the ability to derive solutions over multiple time steps by an LP solver, which assumes the known inflow hydrographs. However, the issue of channel routing has been limited to the use of classical Muskingum approach, where the routing coefficients are constant and calibrated on the basis of a particular historic runoff event, without discussing the effect of using these constant coefficients on alternate regulated flows that can be significantly different from the hydrologic events used during calibration of routing coefficients. Larger basins with multiple reservoirs offer more possibilities for finding better solutions using joint reservoir operation, but they require inclusion of hydrologic channel routing as constraints to optimization. Some recent efforts to address this have been reported by Tahiri et al. (2022), although the routing constraints were associated with iterations imposed within an NFA-based model, which is not an ideal tool to represent the dynamics of channel flow routing since flows along flow links must remain the same, requiring modeling of hydrological transformation of flood ways as a sequence of additions and subtractions to otherwise constant flows in flow links. Other recent attempts to provide comprehensive modeling solution platforms are designed by DetaRes (RTC-Tools 2022) available as a Python package, although the particular details on the implementation of the hydrologic routing schemes are dealt with in a cursory manner requiring more information for independent evaluation.

The above routing model is non-linear, so it typically requires a few numerical iterations where O_t is recalculated and the travel T_s is updated using the interpolated from the table of T_s vs. flows, thus allowing the update of routing coefficients in Equation (4). The process can be repeated a few times until convergence based on a desired tolerance criterion is achieved. The structural form of Equation (4) is identical in form to the well-known Muskingum linear routing form shown in Equation (1), while the coefficients are calculated using a different approach. It should be noted that the SSARR routing coefficients C_i listed in Equation (4) sum up to 1.0, which is also the condition for the Muskingum routing coefficients, and which confirms the conservation of volume. The advantage of the SSARR routing method is that the routing coefficients are a function of flow. As such, they are internally updated as the model simulation moves from low flow season to high flow season. Also, the required travel time vs. flow information for a river reach is much easier to obtain, since it requires the average cross-sectional velocity (which is determined at each flow measurement) and the length of a river reach. For very high flows observed during flood peaks, the hourly timing of the peak is observed for the same hydrograph as it propagates downstream between the two observational points, which provides estimates of travel times during high flows when flow measurements are not conducted due to the associated safety risks. Travel times vs. flow can also be obtained from HEC-RAS modeling by assuming various flow rates, or by obtaining the mean cross-sectional velocities by using the Manning's equation are applying them to the length of the entire river reach.

Note that a typical distance the outflows from the reservoir will travel within 1 day is around 85 km (assuming close to 1 m/s of mean cross-sectional stream velocity), which limits the distance where the storage releases will have effects on downstream flows within a day. Most river basins are longer than 85 km, which implies that storage releases made on day 1 will be felt at downstream locations in subsequent days, depending on the distance. These considerations regarding the travel times are valid both for steady-state model runs as well as for modeling that includes hydrological channel routing. In this context, deriving model releases for individual time steps may or may not meet the downstream operational goals in subsequent time steps, since each time step is an independent model solution. The concept of deriving simultaneous solutions to multiple time steps was introduced to overcome this limitation.

Equation (6) denotes the mass balance at network node n, with A_i representing the incidence matrix coefficient that is typically equal to either −1, 0 or 1, with −1 and 1 denoting the incoming and outgoing flows for network arcs that are associated with a given node n. In multiple time step optimization, the mass balance Equation (6) is applied to each node n in the network and to each time step t over the solution horizon of m time steps. Expression (7) represents the upper bound on flows, which can be either set to constants U_i,t representing for example the maximum storage or canal capacity that should not be exceeded, or it can be defined as a function of flow in some other network component k, such as for example the return flow from irrigation being a linear function of consumptive use or the maximum outflow through an outlet structure such as an orifice or hydro power plant being a piece-wise linearized function of the available storage for a simulated time step.

A complete list of components that can be included in the WEB.BM model is available in previous publications (Ilich 2021). It should be noted that both reservoir and hydrological channel routing are included directly in the optimization process, with only the channel routing constraints requiring iterative calls of the optimization solver, as explained in the ‘Implementation Procedure’ section.

Due to significant reliance on the use of reservoir rule curves, many practitioners consider them to be a useful mechanism for reducing flood damage. However, rule curves were primarily introduced to reduce the chance of premature depletion of storage during dry years. They typically have a fixed shape for each year, and they are unable to offer necessary flexibility for managing floods in large multi-reservoir systems, where each incoming flood has a different spatial distribution and different timing of the peaks, thus requiring a unique response from the system that cannot be obtained from a fixed rule curve shape. The purpose of using optimization is to find the best response within a limited time horizon that was 2 days long in this study, although the time horizons of 1, 3, and 4 days were also tested to evaluate their impact on the quality of model solutions. Hence, the ultimate goal of this study was not to improve the reservoir rule curves, but rather rely on the model to develop short-term operational strategies based on the established allocation priorities and assumed availability of a short-term runoff forecast. The reasons for using LP within the WEB.BM model are listed below:

• LP is the only optimization solver that guarantees finding the global optimum, due to the convexity of it search space that does not allow the existence of local optima;

• There is no need to discretize decision variables into a finite number of possible states as is the case with DP applications, since LP works with decision variables that are double precision real numbers; and

• There is no curse of dimensionality associated with LP. The WEB.BM can solve heavily constrained problems with over 500,000 variables in less than 11 min, while other researchers reported solving problems with close to 5 million variables within hours, as was achieved in a study aimed to set up water management plan for the entire State of California (CALVIN 2023).

It is common wisdom among Operations Research (OR) practitioners that if a problem can be successfully linearized, LP is the best choice for an optimization solver. It is for the above reasons that LP solutions are used as benchmarks to many new experimental solvers in the OR industry. This is unfortunately not the case among water resources practitioners, where respectable benchmark problems and their solutions have yet to be established.

Setting up the pricing vector P_i is easy with early LP-based applications (especially those based on NFA solvers), since the only requirement is to select the P_i values with a relative difference between water use allocated to various stakeholders. For example, giving stakeholders A, B, and C the priority factors of 1, 2, and 3 would yield the same water allocation as would the priority factors 7, 11, and 32, since in both instances the case of maximization (or minimization of deficits) would allocate water first to stakeholder C, followed by stakeholder B and then stakeholder A. Complications started to emerge with the introduction of constraints, as originally discussed by Israel & Lund (1999), and they were further escalated by introduction of full LP solvers and various piece-wise constraints, along with multiple time step solution mode. In spite of this, selection of Pricing vector P_i that will deliver a desirable solution that follows a prescribed operating policy is not a major problem, and it is still considered as a relatively straight forward task for most systems where operational priorities are clearly defined. Although most dams are built for multi-purpose water use, their original design usually had the prime and secondary purpose, i.e., either irrigation was the prime purpose, with hydro power being generated as a byproduct of irrigation releases, or the other way around (considering the year-round operational priorities). Flood events typically last up to one week, and they are usually associated with operation of hydro power plants at full capacity and absence of any shortage of irrigation supply flows, as is the case in this study, where three operational goals have been included in the model in the order in which they are listed:

• (1)

  Flood control;

• (2)

  Water supply for irrigation and industrial demands; and

• (3)

  Hydropower generation.

With major irrigation canals are located downstream of all reservoirs, it is obvious that the only mechanisms that should cause lowering the storage level below the designed full supply level (FSL) are irrigation demands in the low flow season when storage inflows are lower than the downstream demands and pre-flood drawdown during the monsoon season that would minimize the flows bypassing turbines while simultaneously assisting with keeping downstream flows within the full bank capacity. The selected values of pricing vector P_i are provided in the ‘case study’ section. At this point, it should be noted that P_i values have been determined arbitrarily by the user using several trial model runs to ensure that the model operates as expected. In addition to minimizing flood damage, the model was set up to meet the existing water supply demands and maximize generated power, implying a multi-reservoir multi-purpose optimization.

The issue of automated search for the best values of the pricing vector P_i has been addressed by many researchers in recent decades. This was precipitated by new emerging water demands that were not originally formulated when the reservoirs were initially commissioned, such as for example the need to maintain the environmental flow targets, which emerged as an issue in the past few decades. To address this, researchers introduced multi-objective optimization, delivering results in the so-called Pareto-optimal front with possibly hundreds of long-term model solutions that are considered ‘equally optimal’ with respect to two or more management goals which some researchers eventually incorporated into a more general solution procedure such as for example the Evolutionary Multi-Objective Direct Policy Search (EMODPS) proposed by Giuliani et al. (2016). In contrast, reservoir operators seek simple guidelines to determine releases from storage in the short-term, and their allocation priorities in most instances are usually well defined, especially in North America, where water allocation priorities are based on the age of respective water licences – a system that is not likely to change easily due to legal constraints. Quinn et al. (2019) listed several reasons for the existing gap between academia and practitioners, one of which is the lack of ability on the part of reservoir operators to fully understand the selected mathematical utility functions and the black-box tools, due to the lack of understanding the processes involved in their creation. Solution procedures such as EMODPS attempt to guide storage releases based on revising the operating rules. The work presented here is different in the sense that it removes the need to use the ‘upper rule curve’ for any storage reservoirs, but rather dynamically creates operating rules based on the results of optimization for each forecast horizon based on the starting storage levels, runoff forecasts and current level of water demands. The lower rule is not an issue here due to sufficient amount of water to meet the current irrigation demands using the Standard Operating Policy (although this situation may change in future as water demands continue to grow). Lower rule curves are pertinent for managing droughts, and system optimization over the entire operation season is a valuable tool for constructing both the operating zones that correspond to various levels of water rationing policies (Ilich 2023), but this is not applicable to the case study presented here. Hence, the proposed modeling framework mimics the operators' decision that takes into account the existing priorities, starting storage, runoff forecasts and demands, and moves from one day to another by sequentially updating the decision made for a single day based on the available runoff forecast horizon. The WEB.BM model could replicate the same decision-making process for time steps that are shorter than a day, going down to 1 h, but the selection of shorter time steps has a bearing on the length of the downstream river reaches to enable meaningful channel routing calculation. For example, if 1 h is selected as a desired calculation time step, the length of each river reach should be defined such that the average channel flow requires 2 h to traverse the entire river reach, which would typically limit the length of river reach no more than 8 km, thus resulting in creation of numerous river reaches for systems that are more than 400 km in length.

The process of implementing hydrological channel routing into an LP optimization framework requires repetitive calls of the optimization solver until a desired level of convergence is achieved, or until the previously set maximum number of iterations has been reached. Each call of the optimization solver requires that a set of routing coefficients be precalculated based on the most recent channel flows. The routing coefficients are updated be recalculating travel times using the channel flows obtained from the most recent call of the solver, before initiating the next iterative call of the solver. Typically, the maximum number of 5–6 iterations is sufficient to achieve the accuracy of the assumed and recalculated travel times of less than 1%. User have the option to set the tolerance limit (TL) on accuracy and the maximum number of iterations. If the desired accuracy is not achieved within the maximum number of iterations, the model proceeds with the solution from the last iteration and prints out a warning message. This process is presented by using a pseudo code presented below, where an iteration refers to a single call of the optimization solver routine.

For each time horizon solution {

set the tolerance limit TL to a desired threshold; set the maximum number of iterations Imax; set the initial travel times Ts using initial solution flows for a given time of the year;

for all iterations i < Imax {call the optimization solver to get optimal channel flows and reservoir levels;

calculate travel times Ts’ for all channels based on the optimal channel flows;

if[abs(Ts-Ts’) < TL)] // check the compliance with the tolerance limit TL

return to main program with the final solution for a time horizon;

else

set Ts = Ts’ // update the assumed travel times along each channel;

Ci = F(Ts) // recalculate coefficients Ci based on updated the travel times Ts;

i = i + 1; // update the iterations counter;

}

adopt the best optimization solution achieved after the Imax iterations;

}

It should be noted that optimization solver provides flows through the hydro power plants, reservoir spills and diversion canals for the entire system. In essence, these are all reservoir releases (except for the large diversions at Durgapur Barrage), but they are made through different conduits and for different purposes, and as such they are modeled as independent variables.

The WEB.BM model is a web-based application that uses COIN-OR library of public domain solvers Ilich (2021). The model has three solution modes:

• (a)

  Solutions for a single time step forecast (STO-0);

• (b)

  Simultaneous solutions for Multiple Time Step Solutions for n time step forecast (STO-n); and

• (c)

  Simultaneous multiple time step solutions for the entire year (MTO year by year) or for all years.

While options (a) and (b) can be executed with or without the use of hydrological channel routing, option (c) is usually conducted using weekly or longer time steps on a steady-state basis, without the use of hydrological channel routing.

While the solution involves three simultaneous time steps and assumed foreknowledge of flows over a time horizon of 3 consecutive days (Figure 2), the final adopted flows for the first solution are only applicable to day 1, shaded in black (Figure 3). The ending storage for each reservoir on day 1 is used as the starting storage for day 2, and the solution proceeds over the next 3 days. Once this solution is finalized, reservoir operation for day 2 is adopted as final, its storage at the end of day 2 is assumed to be the starting storage for day 3, and the solution proceeds of the next three-day horizon. In this way, the model mimics the decisions of an operator who would take into account a short-term runoff forecast and operational objectives at the start of each day. In addition to the STO-n solution model, for time steps that are multiples of 1 day, the model can be run using two more modes: MTO for year individual year (i.e., all time steps within one year are solved simultaneously) and MTO for all years, which assumes foreknowledge of runoff conditions through the entire period for which input data are available. This can be useful when studying carry over storage requirements from year to year. These two options are typically used for planning studies based on steady-state solutions with weekly or 10-daily time steps. The default model output is shown in the units of flow for channels and net evaporation on reservoirs, elevation (m) for reservoirs, while hydro power output provides both flow through the hydro power plant and the mean power in MW achieved over a simulated time step. The output can optionally be generated in the units of volume, which converts the storage, net evaporation and all channel flows to volumes, thus making it easier to conduct mass balance checks for reservoir nodes. The Google Maps layer shown in Figure 12 is created by the model user with the drawing tools that are internally available in the WEB.BM, which also provides charting tools for both time series and exceedance plots, as well as the statistical analyses of deficits related to the targets for both in-stream and consumptive use targets.

Located in in West Bengal, the Damodar River Basin is one of the oldest basins in India in terms of its history of early water resources developments that date back to the 1950s. There are currently five multi-purpose reservoirs in the Damodar River Basin, built to provide water supply for industry and irrigation, with three installed hydro power plants, but one of main functions assigned to the reservoir operators is protection. Panchet and Maithon are the largest dams in Damodar basin, with total respective installed hydropower plant capacities of 80 and 60 MW. Four out of five reservoirs are currently operated by the Damodar Valley Corporation, and one by the State of Jharkhand. In this paper, the assumption was that all five reservoirs are operated as a system so as to best manage the historic flood events based on the 2-day forecast horizon and the original reservoir design specifications, although some of the reservoirs are currently operated below the originally intended range due to unresolved resettlement and land acquisition issues.

Given a short-term runoff forecast of only 2 days, the optimization program can quickly search through all the available options that also include pre-flood drawdowns. Effective use of the flood storage zones combined with the effects of channel routing and local runoff downstream of Durgapur Barrage can be significant. The principal advantage of LP is that it guarantees finding the global optimum, which is not the case with the heuristic search engines or non-LP solvers.

Based on the established priorities in the Damodar Valley Corporation which is in charge of managing the river basin, hydro power is to be generated as a byproduct of irrigation releases. Had there been no irrigation and no floods, the obvious optimal solution would be to maintain the FSL at all times, in an effort to route all storage inflows through the turbines while maintaining the design operating head. Only one of the following conditions should cause deviations from the FSL:

• (a)

  Model derived drawdown to manage the incoming flows that exceed the hydro power plant capacities and route them all through the turbines by minimizing spills bypassing turbines;

• (b)

  A drawdown caused by irrigation demands being higher than storage inflows; and

• (c)

  A pre-flood drawdown followed by the use of flood storage zone in an effort to maintain the flows in designated downstream river reaches below the full bank capacity.

An obvious optimal model solution is the one that removes flood damage from all historical floods and meets all irrigation water requirements, while maintaining the storage levels as high as possible while minimizing unnecessary spills so as to maximize generated hydro power subject to the above operational constraints. Optimality of model solutions is often hard to establish if arbitrarily selected user-defined utility functions are optimized. When comparing outputs from different models, the relevant optimality criteria should include the number and magnitude of violations of the stated operating goals. Table 1 shows the key reservoir levels used in this study.

The diagram in Figure 13 should be related to the net head which is defined as the difference between the head and tail water elevations, rather than the reservoir water level, but in many instances where tail water elevation is modeled using a constant average for median flows (as is the case in this study), the net head elevation translates to the storage levels of the upstream reservoir. It is obvious that low storage levels cannot produce maximum power, with the typical reduction being in the order of 20–30% of the installed power capacity. The WEB.BM can model this constraint easily using piece-wise linearization of the maximum flow function which can be derived from the net head-flow-efficiency diagrams that are provided by turbine manufacturers. It should be noted that numerous recent publications that detail the use of various heuristic solvers for modeling hydro power typically ignore this constraint.

Flood damage occurs when the flow exceeds 3,000 m³/s at Durgapur Barrage. It is encountered only in the optimization run with a 2-day lead time, and it happens in the following years 1978 (3,092.8 m³/s); 2000 (3,263.5 m³/s); 2006 (3,494.5 m³/s) and 2017 (3,888.7 m³/s). The objective function has a progressive penalty of 9,000 for each unit of flow between 3,000 and 3,200 m³/s, 25,000 for each unit of flow between 3,200 and 3,500 m³/s, and 150,000 for each unit of flow above 3,500 m³/s. Table 2 summarizes the results of four scenarios that differ only by the length of the forecast horizon which was set to 2, 4, 6, and 10 days.

It is also evident from Table 2 that scenarios with 4, 6, and 10 days forecast horizon provide solutions with no flood damage and no deficits to irrigation, in spite of the total annual requirement of 3.2 billion m³. The amount of generated power increases with the increase in the forecasting horizon, until the maximum is found with the required forecasting horizon of 20 days. Increasing the length of the forecasting horizon beyond 20 days does not increase the generated power any more, and therefore represents the best possible solution, as it confirms optimality of all scenarios subject to the stated constraints and the length of the forecasting horizon. These results demonstrate that the combined use of optimization and reliable short-term forecasts can effectively replace the current practice of using the reservoir rule curves during the high flow seasons.

This paper demonstrates that managing reservoirs using optimization can deliver outstanding results provided that accurate runoff forecasts are available. These results indicate that future efforts should focus on improving the runoff forecasting models that should rely on better weather forecasts, remote sensing and possibly the use of artificial intelligence for improved transformation of rainfall to runoff. The synergy of improved forecasts and mathematical optimization for reservoir operation holds out a promise of delivering significant improvements that would save lives and reduce property damage resulting from severe floods. It should also be noted that the WEB.BM can be used both as an operational or a planning model, with a range of time step lengths, and a combination of different multiple objectives that could simultaneously involve maximization of hydro power or minimization of water supply deficits.

Multi-purpose multi-reservoir optimization presented in this paper is easy to understand and verify. It would be useful to compare other similar applications in river basins where all relevant input data could be posted in public repositories, thus enabling cross comparison of results that were obtained with different tools while the input data, constraints and stated objectives remain the same. Obtaining equally good or better results using the EMODPS or other similar modeling approaches compared to multi-purpose multi-reservoir optimization models can help gain traction and interest among practitioners for alternative technologies. Since most modeling approaches claim to provide optimal guidance to reservoir operators for a given set of state variables and short-term runoff forecasts, it would make sense to compare their results. The need for establishing rigorous benchmarks in the water resources modelling sector has been identified in the past (Ilich & Manglik 2022) as a useful guideline to water resources practitioners, especially for water resources systems with moderate complexity.

This project was part of a consultancy agreement funded jointly by the Government of India, Ministry of Jal Shakti and the World Bank under the National Hydrology Project office in New Delhi. The WEB.BM was created by Optimal Solutions Ltd from Calgary, Canada.

The WEB.BM model is accessible free of charge at www.riverbasinmanagement.com requiring only an email address and a password for registration and access, resulting in the creation of a personal workspace within the application for each registered user. Signing on to the platform provides instant access to instructional videos and the User's Manual. Data used in this study belong to the Ganga River Basin, and as such they are restricted from public sharing. However, the WEB.BM User's Manual contains sample data from other basins that can be used to test the capabilities presented in this paper.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.