AbstractThe increasing access to inexpensive sensors, computing power, and more accurate forecasting of storm events provides unique opportunities to shift flood management practices from static approaches to an optimization-based real-time control (RTC) of urban drainage systems. Recent studies have addressed a plethora of strategies for flood control in stormwater reservoirs; however, advanced control theoretic techniques are not yet fully investigated and applied to these systems. In addition, there is an absence of a coupled integrated control model for systems composed of watersheds, reservoirs, and channels for flood mitigation. To this end, we developed a novel nonlinear state-space model of hydrologic and hydrodynamic processes in watersheds, reservoirs, and one-dimensional channels. The model was tested under different types of reservoir control strategies based on real-time measurements (reactive control) and predictions of the future behavior of the system (predictive control) using rainfall forecastings. We applied the modeling approach in a system composed of a single watershed, reservoir, and channel connected in series for the observed rainfall data in San Antonio. Results indicate that for flood mitigation, the predictive control strategy outperforms the reactive controls not only when applied for synthetic design storm events, but also for a continuous simulation. Moreover, the predictive control strategy requires smaller valve operations while still guaranteeing efficient hydrological performance. From the results, we recommend the use of the nonlinear model predictive control strategy to control stormwater systems because of the ability to handle different objective functions, which can be altered according to rainfall forecasting and shift the reservoir operation from flood-based control to strategies focused on increasing detention times, depending on the forecasting.IntroductionFloods are the deadliest natural disaster in the United States and throughout the world (Wing et al. 2020). Estimated global flood damages from 1980 to 2019 exceed USD 750 billion, with a peak in 2012 of nearly USD 70 billion (Our World in Data 2021). Storm events are expected to become more frequent and intense because of climate change, likely increasing not only economic but also social and environmental impacts, posing flood control as one of the greatest challenges for future planning and management of water resources (Gasper et al. 2011). Flood control measures in urban stormwater infrastructures are typically performed by static operations of valves, gates, pumps, and/or tunnels based on predefined heuristic rules. With the advances in real-time control strategies such as the advent of inexpensive sensors, wireless communication, microprocessors, and microcontrollers, opportunities to enhance flood management are evident. Therefore, control theory methods can be applied to enhance water resources management by deploying optimization-based control algorithms. Despite the fact that control theoretic methods have been applied to control combined sewer systems, reservoirs, and drinking water systems (Duchesne et al. 2001; Troutman et al. 2020; Wang et al. 2020), it is a relatively new technique for drainage systems with separated infrastructure for stormwater and sanitary sewage (Wong and Kerkez 2018; Lund et al. 2018).Most drainage infrastructure in major cities was built to operate as static systems. These systems, to convey and/or store large storms (e.g., 100-year storms), typically require relatively large dimensions. However, over the life span, the aging infrastructure, lack of proper maintenance, or the increasing of expected surface runoff (e.g., climate change and urbanization) can decrease the system reliability and hence increase the risk of flooding (Kessler 2011; Zhang et al. 2018). Real-time control (RTC) of the existent systems (e.g., watersheds, reservoirs, channels) can change the flow-storage regime by controlling actuators such as valves and pumps and ultimately restore or increase the level of protection against flooding. However, new stormwater systems designed for RTCs could require smaller surface areas and volumes, potentially leading to an overall cost reduction for the same level of expected performance (Wong and Kerkez 2018; Brasil et al. 2021; Xu et al. 2021).Literature ReviewThe reviewed literature shows several applications of optimization of flow characteristics in a hydraulic structure (i.e., reservoir, channel, pipe, tunnel) to provide multiple benefits and increase the average performance for different water-related problems. An optimized control of the drainage facilities can enhance erosion control (Schmitt et al. 2020), provide stormwater runoff treatment because of higher detention times (Sharior et al. 2019), increase navigability conditions in canals (Horváth et al. 2014), and not only reduce flood downstream locally but also reshape the hydrographs in a desired way set by the developed optimization problem (Wong and Kerkez 2018).Another example of RTC approaches applied to water systems is on the topic of combined sewer overflows, which has an extensive literature (García et al. 2015; Ocampo-Martínez and Puig 2010; Joseph-Duran et al. 2015). However, only a few cities worldwide have applied optimization-based RTC in their drainage systems (Lund et al. 2018). Nonetheless, only a few studies assessed the benefits of RTCs in separated drainage systems. In general, there are two types of control strategies: reactive controls (i.e., based on real-time measurements), and predictive controls (i.e., based on predictions of the future behavior of the system). Reactive controls can be used with heuristic or rule-based approaches and with optimization-based approaches. Predictive controls are typically solved with optimization-based approaches. Recently, Schmitt et al. (2020) assessed the role of heuristic and reactive control rules (i.e., based on real-time measurements) applied to valves in stormwater reservoirs (i.e., actions are made based on measured states), focusing on control of erosive flows. The RTC efficiency was assessed through flow–duration curves, providing exceedance probabilities for any given flow. Although the RTC application provided a significant flow reduction for relatively small flows, its implementation for larger storms (i.e., > 1 year) had slightly less or equal peak flow reduction as the passive control. In several assessed storms, the outflows were larger than the inflows, indicating the RTC increased the likelihood of large flows. Because the controls were guided to increase water quality by increasing detention times (e.g., the control algorithm principle relies on storing the water for 1 or 2 days for any particular event), when sequential storm events occur flood risks increased because the storage capacity was nearly full from previous storms.The research conducted by Wong and Kerkez (2018) presents an optimization methodology for the control of stormwater reservoirs to identify optimal but reactive control based on measured states. A Stormwater Management Model (SWMM) (Rossman 2010) was used to estimate the flow dynamics in a watershed with several reservoirs and links. The model outputs were used in a linear quadratic regulator (LQR) control to decide the valve openings schedule in a set of controlled assets. Although the control schedule given by the LQR reduced peak flows in several ponds, a few had higher outflow peaks in some of the most intense storms assessed, mostly because of the lack of predictability of the future states of the system. Their results illustrate how RTC can increase flood control performance of urban drainage systems, but also indicate that without predictions, its practicality might be limited.In contrast to heuristic and reactive controls, Shishegar et al. (2019) developed an optimization-based approach using a model predictive controller (MPC) for operating valves in a stormwater detention pond. They used a calibrated SWMM model for the watershed and formulated the water balance dynamical problem in the reservoir by linear programming with flows as decision variables. Their approach reshaped hydrographs out of the stormwater reservoir according to the control objective, although several simplifications were adopted: (1) perfect 48-hour rainfall forecasting was used in the prediction horizon, (2) evaporation was neglected, and (3) linear hydraulics released outflows in a linear fashion.Other recent applications of RTC approaches using heuristic controls for water quality enhancement such as the on/off and detention control (Sharior et al. 2019), fuzzy-logic and data-driven algorithms with genetic algorithms (Li 2020), and deep learning (Mullapudi et al. 2020) are found in the literature and address different applications of RTC of urban drainage than our study. From these studies, we categorize RTC for separated urban drainage systems into (1) static and/or optimization-based reactive controls (i.e., control algorithm according to measured or estimated states); and (2) predictive- or optimization-based controls (i.e., control algorithms that consider future states estimation using rainfall forecasting and hydrological models). Several rule-based control (RBC) algorithms do not even require hydrological modeling, although it is mandatory for predictive controllers (Lund et al. 2018).Paper Objectives and ContributionsWe observe from the aforementioned studies a lack of coupled model for the main processes related to floods such as overland flow, infiltration, reservoir routing, and channel routing that allows the use of advanced control theory techniques. Although SWMM or the Gridded Surface Subsurface Soil Analysis (GSSHA) are capable of solving the usually complex shallow-water equations, their use in optimization for large systems can be intractable for real-time control and a simplified plant model of the flow dynamical model is needed (Lund et al. 2018). Moreover, exploring the differences between water quality controls such as detention control and on/off control, further discussed in the next sections, in flood control performance are not yet investigated in the literature. In this paper, we address these issues and explore the trade-offs between reactive and predictive control strategies in regards to flood mitigation. A schematic of the modeled system is shown in Fig. 1, representing the three modeled systems: watersheds discretized in cells, reservoirs receiving outflows from watershed, and channels discretized into subreaches receiving outflows from reservoirs. In this model, we are only interested in floods generated by excess of overland flow. Therefore, subsurface flow are not considered and stored volumes are assumed as overland flow volumes.To this end, we developed a novel state-space representation of the main processes of the water cycle related to urban catchments (i.e., infiltration, overland flow, reservoir storage, channel routing) using cells, reservoirs, and subreaches of channels connected as networked dynamical systems. This state-space model is based on energy, continuity, and momentum equations. We approached the nonlinearities of the flow dynamics by performing successive linearizations in each time step of the model using data from previous time steps as operational points. Although we were able to linearize most of the system’s equations based on the laws of physics (e.g., water balances, energy conservations), we were unable to do it for the rainfall intensity because of the complexity and absence of differentiable models. Therefore, we assumed a known rainfall input time series in the model and we modeled the watershed as a nonlinear dynamical system.In addition to the development of the novel state-space model, we want to assess how varied valve control strategies behave in the system. To that end, we developed an MPC algorithm to enhance the operation of valves in stormwater reservoirs, minimizing a composed cost function related to flood performance. Moreover, we compared the efficiency of reactive controls (i.e., controls based on real-time measurements of the states) with MPC solved with a gradient-based method (interior-point).The fundamental contributions of this paper are described as follows: •We present an overall mathematical representation including the flow dynamics in watersheds using the nonlinear reservoir and the Green-Ampt infiltration model (Green and Ampt 1911), coupled with reservoir routing and one-dimensional (1D) channel dynamics in a nonlinear state-space representation. This allows optimization and real-time control of urban drainage systems without requiring extensive software packages (i.e., only MATLAB is required).•We derive the nonlinear dynamics of watersheds and linearize the dynamics of reservois and channels. Because no continuous and differentiable model is currently available for rainfall intensity, we assume a known time series of rainfall as piecewise continuous input data.•We provide a comprehensive analysis of reactive controls (i.e., some rule-based and other optimally controlled) for flood mitigation and tested the efficiency of water quality rule-based algorithms presented in Sharior et al. (2019) for flood mitigation.•We develop and apply a servo control algorithm (Young and Willems 1972) used in the discrete linear quadratic integrator (DLQI) reactive control. This is a new application for this control technology in urban drainage systems. This algorithm allows tracking a specific state and can be used in reservoirs with minimum specified water surface depths (e.g., wetlands or retention ponds).•We evaluate and discuss the performance of reactive RBCs (passive, on/off, detention control) and reactive optimization-based controls [discrete linear quadratic regulator (DLQR) and discrete linear quadratic integrator] compared with the predictive control strategy. Moreover, we explore the caveats of consecutive design storms and also compare the flood performance of reactive and predictive algorithms in a continuous simulation, providing a methodology to assess the efficiency of control algorithms for flood mitigation.The remainder of the paper is organized as follows: “Mathematical Model Development” develops the state-space model for watersheds, reservoirs, and 1D channel flow-routing dynamics. Next, “Reactive and Predictive Control Strategies” describes the tradional reactive controls applied to control of stormwater reservoirs. Moreover, in this section we develop a novel nonlinear MPC optmization problem focusing on flood mitigation in channels and reservoirs. Following, in “Mathematical Model Application” we present a case study to test the controls presented in the aforementioned sections, including two scenarios: two consecutive design storms of 25 years, 12 hours and 10 years, 12 hours, respectively, and a continuous simulation scenario from April 23 to July 23, 2021 in San Antonio. “Results and Discussion” shows the results and discussion of the model application and “Conclusions, Limitations, and Future Works” the conclusions, limitations, and future works. The notation for this paper is introduced next.Italicized, boldface upper and lowercase characters represent matrices and column vectors: a is a scalar, a is a vector, and A is a matrix. Matrix In denotes an identity square matrix of dimension n by n, whereas 0m×n and 1m×n denote a 0 and 1 matrix with size m by n, respectively. The notations R and R++ denote the set of real and positive real numbers. Similarly, N and N++ denote the set of natural and positive natural numbers. The notations Rn and Rm×n denote a column vector with n elements and an m-by-n matrix in R. The element-wise product or Hadamard product is defined as x○y≔[x1y1,x2y2,…,xnyn]T multiplications. Similarly, the element-wise division or Hadamard division is defined as x⊘y≔[x1/y1,x2/yn,…,xn/yn]T. The element-wise p power of a matrix A, (A°p), with A∈Rm×n and p∈R is given by ai,jp for i∈N++ and j∈N++ The number of elements in a set A∪B is n(A∪B)=n(A)+n(B)−n(A∪B). A normally distributed random number with average μ and variance σ2 is notated by N(μ,σ2) Given a vector x∈Rn, the notation x(i:j) with i and j∈N++ represents a cut in x from the ith to jth entries.Mathematical Model DevelopmentThe stormwater flow dynamical problem solves physics-based governing equations in each watershed, reservoir, and channel using physically based input data. Because mass balance equations are solved, we postulate the stormwater flow dynamical system as a nonlinear difference-algebraic equation (DAE) state-space model. All variables used in this paper are summarized in Table S1.The model requires matrices and vectors to represent the dynamical and algebraic parts of the simulated hydrological systems, as presented in Eq. (1). Specifically, matrix E enables the representation of the dynamics and algebraic constraints of the coupled watershed-reservoir-channel model in a single state-space model. Moreover, linear time-varying parts are represented by matrices A(k),B(k), while C represents time-invariant output matrix because the sensors are assumed to have fixed geographic placements. In addition, offsets, nonlinearities from the watershed model, operational points, rainfall intensity in each cell, outflows and inflows connectivity from watersheds to reservoirs and from reservoirs to channels, and integrator reference setpoints are given by ψ(·). The mathematical development of these matrices and vectors is detailed in the following sections. In this paper, we develop the DAE state-space model for a system composed of a single watershed, reservoir, and channel. We collect a vector of water surface depths in cells, reservoirs, and channels, accumulated infiltration depths in each cell, and outflows from catchments and reservoirs as the state vector, such that x(k)=[hefw(k),fdw(k),qoutw(k),hr(k),qoutr(k),hc(k)]T. We also assume a control vector given by u(k)=ur(k). The model parameters, states, outputs, and sources of uncertainty are presented in Table 1. The state-space representation can be written as (1a) Ex(k+1)=A(k)x(k)+B(k)u(k)+ψ(k,x(k),x*r,u*r)(1b) where E = singular matrix with some zero rows representing the algebraic constraints of flow equations; A(k)∈Rn×n = state or system matrix; x(k)∈Rn = state vector; B(k)∈Rn×m = input matrix; u(k)∈Rm = input vector; ψ(k,x*r,u*r) = disturbance vector; x*r and u*r = operational points that are updated each time step; C∈Rp×m = output matrix; y(k) = output function; and k = time step index. The vectors hefw,hr, and hc are water surface depths in each system, fd is the accumulated infiltration depth in the cells of the watershed, and qoutw and qoutr are the watershed outflow and reservoir outflow, respectively. In the following sections, we define each system and their governing equations as well the linerizations.Table 1. Systems, parameters, states, outputs, and uncertainty features of the modelTable 1. Systems, parameters, states, outputs, and uncertainty features of the modelSystemParametersStates (symbols)Outputs (symbols)UncertaintyCellsInfiltration parametersWater surface depth in each cell (hef∈Rq) (mm), infiltrated depths (fd∈Rq) (mm), and outflow (qoutw) (m3/s)Outflow in each cell and in the outlet (qout∈Rq) (mm·h−1)Manning’s coefficient and rainfall spatial distributionSurface roughnessInitial abstractionDigital elevation modelLULCReservoirsStage-discharge relationshipsWater surface depth (hr∈Rnr) (m), and outflow qoutr (m3/s)Maximum water surface depth (max(hr)∈R) (m)Water surface depth measurement noiseArea-volume function porosityChannelsStage-discharge relationshipsWater surface depth (hc∈Rnc) in each subreach of the channel (m)Maximum water surface depth (max(hc)∈R) (m)Manning’s coefficientManning’s coefficientHydraulic radiusGridded bathymetryFor cases with more watersheds, reservoirs, and channels, the state vector can be augmented to include the new systems, concatenating each individual state (e.g., hefw(k)=[hef,1w,…,hef,sw]T) such that x(k)=[hefw(k),fdw(k),qoutw(k),hr(k),qoutr(k)]T, u(k)=[u1r(k),…,usr(k)]T, where s is the number of systems composed by watersheds, reservoirs, and channels. In this paper, we present the mathematical formulation for a watershed-reservoir-channel system.Watershed Overland Flow ModelingIn this section we derive the mathematical formulation to estimate overland flow in watersheds using a fully distributed hydrologic model. Excess of infiltration (Hortonian flow) and/or saturation (Dunnian flow) generates overland flow in storage cells (Maxwell et al. 2014). In urban environments with high impervious areas, Hortonian flows govern the overland flow generation and occur when infiltration capacity is smaller than the inflow rate (i.e., net precipitation, inflow from neighbor cells, ponding depth).The infiltration losses are estimated using the Green and Ampt infiltration model. This model is physically based and derived from simplifications of the Richards’ equation (Richards 1931; Green and Ampt 1911). All parameters of the model can be estimated in laboratory tests. However, substantial studies are available in the literature providing good parameter estimates according to the soil characterization (Green and Ampt 1911; Rossman 2010). Furthermore, Green-Ampt soil properties are typically available in spatially distributed geographical information system (GIS) databases in many parts of the world, which facilitates the application of this model.For each cell of a predefined grid domain, the (1) saturated hydraulic conductivity, (2) suction head pressure (capillarity), (3) initial and saturated soil moisture, and (4) initial infiltrated water content are required. Extractions of ASCII files from delineated watershed rasters can be used to define the matrices representing the digital elevation model (DEM) and the imperiousness map. More details of these files and the model construction are found in Gomes et al. (2021a, b). Manning’s equation is used to relate water surface depth to flow for the cells (Akan 1993; Te Chow 2010). To perform the calculations, input data such as (1) Manning’s coefficient; and (2) central elevation of each cell of the grid are required.One-Dimensional Vertical Infiltration ModelThe Green-Ampt model is applied to each grid cell to estimate the infiltration capacity at any given time and is used to estimate the available depth to be routed to downstream cells. The infiltration capacity is estimated as (2) ci,j(t)=ksati,j[1+(ζi,j+hefi,j(t))(θsi,j−θii,j)fdi,j(t)]where the subindexes i and j = cell position in the grid; c(t) = infiltration capacity (mm·h−1); ksat = saturated hydraulic conductivity (mm·h−1); ζ = suction head pressure (mm); hef(t) = water depth in the cell (mm); Δθ=θs−θi is the effective soil moisture; and fd(t) = time-varying accumulated infiltration (mm).The infiltration model is a nonlinear time-varying function of the accumulated infiltrated volume Iii,j(t) and is dynamically computed in an explicit first-order finite-difference discretization, written as (3) fdi,j(t+Δt)=fdi,j(t)+fi,j(t)Δt=fdi,j(t)+[min(ci,j(t),qini,j(t)+ipi,j(t)−eTRi,j(t))]⏞fi,j(t)Δtwhere ip = rainfall intensity (mm·h−1); eTR = real evapotranspiration intensity (mm·h−1); Δt = model time step; fdi,j(t) = accumulated infiltration depth (mm); qini,j(t) = inflow discharge rate (mm·h−1); and hefi,j = runoff water depth (mm) in cell i, j.In the previous equation, we model the soil drying only by assuming a flux of evapotranspiration occurring at pervious surfaces such that in drying periods the soil storage depth is decreased. Moreover, we limit the accumulated infiltrated soil depth fd to a minimum threshold typically assumed as 5 mm.Vertical Water BalanceHydrological processes of infiltration, precipitation, surface runoff, and evaporation occur simultaneously. However, an analytical solution of the continuous functions of these inputs into the overall water balance is typically not available for real case scenarios, especially because of the rainfall. With proper stable time-step resolution, an alternative to solve the overland flow dynamics is to derive an explicit system of equations from the time derivative of the storage variation in each cell, expressed as follows (Rossman 2010): (4a) dsi,j(t)dt=[qini,j(t)+ipi,j(t)−eTRi,j(t)−qouti,j(t)−fi,j(t)]ωw(4b) dhefi,j(t)dt=1ωwdsi,j(t)dtwhere si,j = ponded water storage (m3); ωw=ΔxwΔyw is the cell area (m2); Δxw and Δyw = cell resolution in x and y (m); qouti,j(t) = outflow discharge rate (mm/h); and fi,j(t) = infiltration rate (mm/h).Generalizing Eq. (4b) for a vector notation concatenating the number of rows and columns into a vector of dimension q, changing the units for water surface depth instead of storage, we can derive the following expression: (5) hef(k+1)=hef(k)+Δt600(qin(k)+ip(k)−eTR(k)−qout(k)−f(k))where k = time-step index, hef, qin, ip, eTR, qout, and f∈Rq.Outflow DischargeThe outflow discharge is simplified to be a function of the steepest 8-direction slope, cell roughness, and water surface depth in a kinematic-wave shallow-water simplification approach. This approximation is implemented because the goal of the watershed model is to determine flows and not high-resolution surface water flood depths. According to the gridded elevation, first we define a flow direction matrix similar as presented in Fig. 2(b). This matrix is determined calculating the steepest topographic slope of each cell assuming eight boundary cells (i.e., Moore neighborhood grid). According to the steepest slope direction, a number is assigned to each cell. This matrix defines the boundary conditions among each cell in the grid.Using Manning’s equation (Te Chow 2010) and assuming the energy slope as the bottom slope, we can estimate the outflow discharge for a specific cell as follows in a matrix notation: (6) qout(k)=(kfΔxw+Δyw2)s0°1/2⊘n⏞λ∘(max(hef(k)−h0),0)°5/3where kf = conversion factor equal to 1×10−5; s0 = steepest topographic slope (m/m); n = Manning’s roughness coefficient (sm−1/3); λ lumps the hydraulic properties of the cells into a single vector; and s0, n, h0, qout, and λ∈Rq. All the other equations are in the international system of units.To guarantee continuity in the hydrological model, a preprocessing in the digital elevation model is performed. Natural sinks are filled to ensure all cells have an outlet slope. Moreover, the outlet boundary condition of the watershed is modeled assuming a normal flow hydraulic condition (Kollet and Maxwell 2006).Inflow DischargeThe inflow discharge is a function of the flow direction matrix and outflow discharge. Defining a matrix Bdw∈Rq×q to represent the direction relationship among the cells in the form of a sparce matrix filled with ones (i.e., outflows becomes inflow for the downstream cell) and zeros [i.e., no flow connection between cells, see Fig. 2(c)], we can write (7) where Bq∈Rq×q is the direction Boolean matrix containing the relationship between each cell.Substituting Eqs. (7) and (6) into Eq. (5), tracking the accumulated infiltration depth (fd), and the watershed outflow (qoutw), the watershed subsystem from Eq. (1) can be written in a state-space representation, given by (8) [I2q000]⏞Ew[hef(k+1)fd(k+1)qoutw(k+1)]⏞xw(k+1)=I2q+1⏞Aw(k)xw(k)+Δt[1/600((Bdw−Iq)∘λ∘max(hef(k)−h0,0)°5/3+ip′(k)−f(k))f(k)−0.277×10−6Δt∥λi0,j0max(hefi0,j0(k)−h0,0)°5/3ωw∥]⏞ψw(k,xw(k)where xw(k)=[hefw(k),fdw(k),qoutw(k)]T; i′(k)=ip(k)−eTR(k); f(k) is modeled in Eq. (3), i0 and j0 = indexes of the outlet cells; and λi0,j0 concatenates λ for outlet cells.Reservoir DynamicsIn this section, the reservoir routing dynamics are described and we provide a fully linearizable model to account for valve control in stormwater reservoirs with hydraulic devices as orifices and spillways.Orifice ModelingThe control signal u(k) represents the percentage of the orifice area that allows flow to be routed to downstream channels. Therefore, applying the energy equation in the reservoir (Te Chow 2010) and including u(k) in the effective orifice area, we can derive the controlled orifice equation, such that (9) qo(hr(k),u(k))=ur(k)cd,oao2g(max(hr(k)−(ho+hm),0)=ur(k)koh^r(k)where qo = orifice discharge; cd,o = orifice discharge coefficient; ao = orifice area; g = gravity acceleration; h^r(hr(k))=max(hr(k)−(ho+hm),0) is the effective water depth at the orifice; u = control input representing the valve opening between 0 and 1; ho = bottom elevation of the orifice; and hm = minimum water surface depth to begin the outflow (i.e., typically 20% of the hydraulic diameter of the outlet) (Te Chow 2010).Spillway ModelingThe spillway is also assumed to be discharging at the atmospheric pressure. The Francis spillway equation is typically used for detention reservoirs and can be modeled as follows: (10) qs(hr(k))=cd,slef(hr(k)−p)3/2=ks(hr(k)−p)3/2where qs = spillway discharge; cd,s = spillway discharge coefficient; and lef = effective length of the spillway (Te Chow 2010).Reservoir OutflowThe outflow in a reservoir is a function of the water surface depth hr(k) and is described by the energy conservation applied into the orifice and spillway, and thus has two governing equations. The first case is where hr(k)≤p and is given by Eq. (9) (Te Chow 2010). When the water level reaches the spillway level, the reservoir outflow (qoutr) is the sum of the orifice and spillway flow. Therefore, the reservoir outflow function be derived as follows: (11) qoutr(hr(k),ur(k))={ur(k)koh^r(k)if hr(k)≤p,elseur(k)koh^r(k)+ks(hr(k)−p)3/2We compute the Jacobian matrix of qoutr with respect to hr and u to obtain a linearized flow equation neglecting the high-order terms of the Taylor’s series, resulting in the following equations: (12) ∂qoutr(hr(k),ur(k))∂h≈{ur(k)ko2h^r(k)+3ks[max((hr(k)−p)1/2,0)]2}=α(hr(k),ur(k))(13) ∂qoutr(hr(k),ur(k))∂u=koh^r(k)=β(hr(k),ur(k))Therefore, a linearized model for the outflow in terms of the stored water surface depth and valve opening is given as (14) qoutr(hr(k),ur(k))=qoutr(h*r)⏞γ(k)+α|h=h*r,u=u*r⏞α˜(k)(hr(k)−h*r)+β|h=h*r,u=u*r⏞β˜(k)(ur(k)−u*r)where γ = offset; α˜ = linear coefficient with respect to hr; β˜ = linear coefficient in terms of u; and u*r and h*r = operation points given by the states and controls of the previous time step, such that u*r(k)=ur(k−1) and h*r(k)=hr(k−1).Reservoir Water BalanceThe temporal evolution of storage in a reservoir depends on the inflow, precipitation, evaporation, water surface area, and stage-discharge function of the outlet hydraulic devices. Applying evaporation and precipitation in the reservoir surface area, we can derive an expression for the water storage dynamics, given by (15) ∂sr(hr(k),ur(k))∂t=qoutw(k)+(i(k)−ev(k))ωr(hr(k))⏞qinr(k,hr(k))−qoutr(hr(k),ur(k))where sr = stored volume of stormwater runoff in the reservoir; ωr(hr(k)) = reservoir surface area in terms of hr(k); ev = evaporation in the reservoir surface area; qoutw = inflow from upstream catchment; qinr = total inflow; and i = rainfall intensity. Assuming an average porosity η representing the stage-storage relationship in the reservoir (e.g., for free surface reservoirs, the porosity equals 1), the water surface depth dynamics can be derived as (16) ∂h(k,hr(k),ur(k))∂t=1ωr(hr(k))η[qinr(hr(k))−qoutr(hr(k),ur(k))]The storage dynamics in a reservoir can be very slow depending on the area of the reservoir, which might contribute for low degrees of controllability for reactive controls, even in events with high inflows. Substituting the linearized reservoir outflow, Eq. (14), into the water balance equation, it follows that (17) ∂h(hr(k),ur(k))∂t=1ωr(hr(k))η⏞μ(hr(k))[qinr(k,hr(k))−α˜(k)(hr(k)−h*r)−β˜(k)(ur(k)−u*r)−γ(k)]Assuming an approximated finite-difference scheme by the forward Euler method applied in the water surface depth partial derivative equation, we obtain (18) ∂h(hr(k),ur(k))∂t≈hr(k+1)−hr(k)ΔtGeneralizing the reservoir dynamics for more than one reservoir per watershed (nr>1 and ur(k)∈Rnr), substituting Eq. (17) into Eq. (18), and expanding for a matrix notation, the water depth dynamics in reservoirs as: (19) hr(k+1)=(Inr−diag(Δtα˜(k)∘μ(hr(k))))⏞A˜r(k)hr(k)+(−diag(Δtβ˜(k)∘μ(hr(k))))⏞B˜r(k)ur(k)+Δtμ(hr(k))∘(α˜(k)∘h*r+β˜(k)∘u*r−γ(k)+qinr(k))⏞ψ˜r(k,u*r,x*r)where hr(k),α˜(k),β˜(k), and ψ˜r(k,u*r,x*r)∈Rnr; nr = number of reservoirs; A˜(k) and B˜r(k)∈Rnr×nr; and qinr(k) captures the watershed-reservoir outflow/inflow connection.Moreover, generalizing the method for s systems with nr reservoirs per watershed, we can define z=nrs resulting in the reservoir state-space nonlinear dynamics tracking the reservoir outflows such that (20) [IzOz×zOz×zOz×z]⏞Esr[ysr(k+1)φout,sr(k+1)]⏞xsr(k+1)=[A˜1r(k)O⋯OOOA˜2r(k)⋯OO⋮⋮⋱⋮⋮OO⋯A˜sr(k)OOO⋯O−Iz]⏞Asr(k)[h1r(k)h2r(k)⋮hsr(k)φout,sr(k)]⏞xsr(k)+[B˜1r(k)O⋯OOB˜2r(k)⋯O⋮⋮⋱⋮OO⋯B˜sr(k)OOOO]⏞Bsr(k)[u1r(k)u2r(k)⋮usr(k)]⏞σsr(k)+[ψ˜1r(k)ψ˜2r(k)⋮ψ˜sr(k)σsr(k)∘ko*∘(h^sr(k))°1/2+ks*∘(max(ysr−p,0))°3/2]⏞ψsr(k)where ko* and ks* collect ko and ks for all reservoirs in all systems, respectively. Similarly, p collects spillway elevations. Matrices Esr and Asr∈R2z×2z, matrix Bsr∈R2z×z, while xsr∈R2z, ψsr∈R2z, and σsr∈Rz. Vectors φout,sr,ysr, and σsr are defined as follows: (21a) φout,sr(k)=[qout,1r,1(k),…,qout,1r,nr(k)⏞qout,1r(k),…,qout,sr,1(k),…,qout,sr,nr(k)⏞qout,sr(k)]T(21b) ysr(k)=[h1r,1(k),…,h1r,nr(k)⏞h1r(k),…,hsr,1(k),…,hsr,nr(k)⏞hsr(k)]T(21c) σsr(k)=[u1r,1(k),…,u1r,nr(k),⏞u1r(k)⋯,usr,1(k),…,usr,nr(k)⏞usr(k)]TOne-Dimensional Channel DynamicsThe flow dynamics in rivers or channels can be modeled using hydrologic or hydraulic routing modeling approaches. The first is governed by the water balances in the inlet and outlet sections of a subreach. The continuity equation and storage-outflow relationships are used to estimate the outflow in the last subreach (McCarthy 1938; Cunge 1969). However, this approach is suitable only for estimating flows. Moreover, hydrologic routing is not flexible enough to explain more detailed phenomena such as backwater effects from downstream reservoirs or flood waves in channels with very flat slopes (Cunge 1969). Another issue is the time-related coefficients associated with the hydrologic routing equations. Typically, to ensure the model’s stability, coefficients are in the order of several hours and days. This poses as a drawback for real-time monitoring of urban channels (Kumar et al. 2011). We propose a diffusive-wave simplification in the Saint-Venant equations (SVEs) to represent the flow dynamics in the 1D channels.Recently, an application of the full one-dimensional Saint-Venant equations allowed the state-space representation in channels (Bartos and Kerkez 2021). The authors developed a backward Euler implicit scheme solving continuity and momentum equations via a sparce matrix system of equations. While model stability is theoretically increased with this implicit numerical scheme, to develop a full state-space representation of watersheds, reservoirs, and channels, a implicit derivation for the other systems would be required. Therefore, we solve the 1D flow dynamics using an Euler explicit numerical scheme as used in Eqs. (8) and (20).Channel Water Depth DynamicsTo represent the channel 1D dynamics in space and time, we develop an explicit diffusive-wave simplification in SVE, assuming the friction slope from Manning’s steady-flow equation. Fig. 5(c) shows a scheme of a 1D channel. The outflow in each segment is calculated through Manning’s equation as follows (Panday and Huyakorn 2004): (22) qoutc(hic(k))=aic(k)(rh,i(k))2/3ni(∂hic(k)∂y)1/2where aic = wetted area; rh,i = hydraulic radius; hic = channel water surface depth at cross section i; and y = channel’s longitudinal direction. The previous expression is expanded in a vector format given by (23) qoutc(hc(k))=1nc×1⊘nc∘ac(hc(k))∘(rhc(hc(k))°2/3○(∂hc(k)∂y)°1/2where nc, ac, and rhc = Manning’s coefficient, cross-section area function, and hydraulic radius for each subreach.A channel is discretized into subreaches with appropriate spatial resolution to be suitable for the time step of the model (Chang and Wang 2002). Therefore, the model simulates steady nonuniform flow assuming initial boundary conditions from qoutr(k) and from the outlet normal friction slope (Panday and Huyakorn 2004; Maxwell et al. 2014). The conservation of mass and momentum equations are given by (24a) ∂qic(k)∂y=−∂aic(k)∂t(24b) ∂hic(k)∂y=s0c−(qic(k)nicai(k)rh,i(k)2/3)2⏞sf,i(k)where y = longitudinal channel dimension; qic = net flow within cross sections; s0c = bottom slope; nic = Manning’s roughness coefficient; and sf,i = friction slope.Channel Water Surface Depth DynamicsWe can expand Eq. (24a), resulting in a vectorized channel water surface depth mass balance as (25) hc(k+1)=hc(k)+Δt1nc×1⊘(Δx○Δy)(qinc(k)−qoutc(k))⏞Δq(hc(k))where hc = water surface depth in each subreach of the channel; and the outflow qoutc is given by the Manning’s equation that depends on known functions of cross-section area and hydraulic radius for each reach of the channel and length (Δy) and width (Δx). For the first subreach, qinc is equal to the reservoir outflow (qoutr). Flows qin,ic(k) and qout,ic(k) are equal to qi−1c and qic, respectively.To explicitly solve Eq. (25), we need to derive an expression for Δqc(hc(k)). To this end, we apply the energy conservation within two consecutive cross sections neglecting the velocity head to develop a vector representation for the friction slope assuming normal depths, such that for a particular cell (i) in a 1D channel domain, we can write (26) sf,i(k)=1Δyi(el,i+hic(k)−el,i+1−hi+1c(k))where el = bottom elevation of the subreach segment; and sf,i = friction slope in cell i at step k.Expanding the previous equation in a vector format, we can derive the momentum equation from Eq. (24b) for all subreaches as a linear combination of hc(k) [Eq. (23)], such that (27) ∂hc(k)∂y=Aslopehc(k)+bslopewhere (∂hc)/(∂y)∈Rnc is a vector representing the water slopes in each reach of the channel; and Aslope∈Rnc×nc and bslope∈Rnc are derived from Eq. (26) assuming 1D connections between cells, computing the water surface slopes, and including the outlet slope boundary condition.One-Dimensional Flow in Open ChannelsSimilarly, to the overland flow model for the watershed model, we can compute (qinc(k)−qoutc(k)) as a function of a flow direction matrix and the flow through Manning’s equation. Therefore, the net flow Δqc(hc(k))=qinc(k)−qoutc(k) for all subreaches can be given as (28) Δqc(hc(k))=Bdcqoutc(hc(k))+w(k)where Bdc = topology matrix linking each subreach segment with the previous one; and w(k) = zero vector where only the first entry is equal to the reservoir outflow (qoutr(k)).However, to estimate the 1D flow propagation in channels, we do not assume a simplified hydraulic radius and define a constant λ because of the relatively small width of the subreaches when compared to the gridded cells from the watershed system (Liu and Singh 2004). To estimate λ, functions that describe the cross-section area (ahc(k)) and hydraulic radius (rhc(k)) are required and can be derived in terms of topographic properties of the channel. For a rectangular channel (i.e., condition of the case study), ahc(hc(k))=Δx∘hc(k) and rhc(hc(k))=ahc(hc(k))⊘(Δx+2hc(k)).Linearized Channel DynamicsWe can substitute Eq. (28) into Eq. (25), resulting in a matrix expression for the channel nonlinear dynamics given by (29) hc(k+1)=hc(k)+Δt1nc×1⊘(Δx∘Δy)∘[Bdc1nc×1⊘nc∘ahc(hc(k))∘rhc(hc(k))°2/3∘(Aslopehc(k)+bslope)°1/2)+w(k)]Finally, we can obtain a linearized expression around an operation point h0c, applying the Jacobian in the previous equation, such that (30) Inc⏞Echc(k+1)=(Inc+∇hc(h0c))⏞Ac(k)hc(k)where ∇ = gradient of Eq. (29) and can be obtained knowing the cross-section area and hydraulic radius functions of each subreach and can be computationally derived with symbolic functions in MATLAB, for instance. The gradient and the operational point from Eq. (30) are refreshed each model time step.State-Space RepresentationThe derivation of the state-space model for a single watershed, reservoir, and channel is performed by merging the dynamical state-space representation for each subsystem from Eqs. (8), (20), and (30), and applying Eq. (1), the watershed-reservoir-channel dynamical system can be defined as (31) [EwOOOErOOOEc]⏞E[xw(k+1)xr(k+1)xc(k+1)]⏞x(k+1)=[Aw(k)O(2q+1)×1O(2q+1)×ncO1×(2q+1)Ar(k)O1×ncOnc×(2q+1)Onc×1Ac(k)]⏞A(k)[xw(k)xr(k)xc(k)]⏞x(k)+Br(k)⏞B(k)u(k)+[ψw(k,xw(k))ψr(k,x*r,u*r)Onc×1]⏞ψ(k)Similarly, the derivation of the state-space dynamics for the general case with multiple watersheds, more than one reservoir per watershed, and more than one channel per reservoir can be done using Eqs. (20) and (8). The watershed is an autonomous system, whereas the reservoir is governed by the control law and the channel depends on the outflow from the reservoir, although it is not directly controlled.Reactive and Predictive Control StrategiesThe control strategies tested in this paper are focused on reducing flood effects. Here in this paper, we define it as a composite function accounting for local floods in reservoirs and channels while minimizing rapid changes in valve operation. We can generally write the flood performance as a function defined for (ns) systems accounting for flooding (e.g., reservoirs and channels) with different numbers of components for each system that impacts the flood performance (c(i)) (e.g., flood performance for reservoirs is associated with water levels and control deviations, whereas for channels it is only associated with water levels). Therefore, we can generally write flood performance as linear combinations between weights and functions for each component of each system associated with floods, such that (32) Flood Performance=∑is=1ns∑jc=1c(i)αis,jcfis,j(X)where X = state matrix concatenating x(k) for all simulation time.The aforementioned equation is defined in detail in the objective function used to describe the flood mitigation performance in the section “Model Predictive Control.” This particular function operates in all modeled states returning a real number. Although the flood performance was defined here, only the model predictive control strategy is able to control all parts of the system (e.g., channels and reservoirs) and act directly to maximize it because the reactive controls are developed only using the reservoir dynamics. This section is organized as follows: first we describe the reactive control strategies tested in this paper, and following that we introduce a nonlinear model predictive control algorithm to control valve operations in stormwater reservoirs.Discrete Linear Quadratic RegulatorGiven a linearized time-invariant state-space representation of the reservoir dynamics, the linear quadratic regulator strategy finds an optimal control signal for time step k minimizing a quadratic cost function, solved as an unconstrained convex optimization problem given by (33) minu  J(u)=∑n=1∞[(hr(k))TQhr(k)+(ur(k))TRur(k)]where Q and R = tuned symmetric positive semidefinite matrices representing the set points, cost of the states, cost of controls, and cross-term costs.This particular control strategy only considers the reservoir dynamics. Therefore, the coupled flood performance of this case only focuses on reservoir water depth and control energy as shown in the aforementioned equation. During the period where the dynamics are unchanged (e.g., steady inflows from the watersheds), the performance of the LQR control is optimal (Wong and Kerkez 2018). The solution of this equation given by the discrete algebraic Riccati equation (DARE) (Pappas et al. 1980) consists of finding a P(k) matrix (cost-to-go function) and a closed-loop feedback gain K(k), such that (34) K(k)=(R+(Br(k))TP(k)Br(k))−1(Br(k))TP(k)Ar(k)where the aforementioned matrices P(k) and K(k) can be solved using the DLQR function in MATLAB.A schematic detail of the closed-loop system for the DLQR is presented in Fig. 3(a). In this control algorithm, we are only using the reservoir plant in the dynamics. To design the state feedback gain matrix K, we can tune Q=CTC and R=Inr as a first attempt. Basically, the DLQR control tends to find an optimal control for time (k+1) that ultimately stabilizes the reservoir dynamics or, in other words, tends to release the water considering the costs of rapid changes in the valves. To ensure physical constraints in the inputs (Wong and Kerkez 2018), we limit the feedback gain to (35a) (35b) 0≤uir(k)≤1,i={1,2,…,nr}Discrete Linear Quadratic IntegratorThe linear quadratic integrator works similarly to the DLQR control; however, an augmented dynamical system to account for a reference tracking set point (Young and Willems 1972) is developed. The system dynamics are augmented to include the temporal evolution of the integral of the error between the reference and the output. The set-point reference can be a constant value or a time-varying function representing the goal of the control technique. Assuming the error between states and outputs as e˙(k)=r(k)−y(k) and discretizing it by an explicit forward Euler method in terms of the error integral (e(k)), the augmented dynamical system can be given as (36) [hr(k+1)e(k+1)]=[Ar(k)Onr×nr−CrInr][hr(k)e(k)]+[Br(k)Onr×nr]ur(k)+ψ˜r(k,x*r,u*r)+Δt[Onr×11nr×1]r(k)where r(k) = set point or reference goal (e.g., minimum water surface depth to maintain natural ecosystems in the reservoir); and Cr=Inr.The closed-loop system for the DLQI is presented in Fig. 3(b). A known input is given by qoutw(k), which changes the dynamics over time. The solution of the DARE is performed with the coefficients of matrix Q [Eq. (36) and Fig. 3] representative of the deviations in states and outputs and from matrix R the control signal of 1×100, 1.5×103, and 1×102, respectively. The number of integrators is assumed to be equal to the number of outputs; hence, all nodes are considered observable. After determining the optimal feedback gain, to ensure physical constraints, we also clip the control signal following Eq. (35) as shown in Figs. 3(b and c). The control input can be separated into (37) ur(k)=−[kske][hr(k)e(k)]where ks = system feedback gain; and ke = servo feedback gain.The ultimate objectives of DLQR and DLQI are presented in Fig. 3(c). While DLQR tends to find a control schedule steering the output function to zero [e.g., rapidly release stored water and decrease hr(k) considering control energy], DLQI has a similar approach but for a reference water level r(k)DLQI (i.e., a given expected reservoir water level).Rule-Based ControlsTypically, stormwater systems are passive or controlled by rule-based reactive and local controls, sometimes even manually by operators (Schmitt et al. 2020; Shishegar et al. 2019; García et al. 2015). Although the control rule of these methods seems logical and easily applicable for simple storm events, they can fail to control more complex cases such as fast consecutive storms by not predicting the future behavior of the system (Lund et al. 2018; Sharior et al. 2019). Some of the most common rules for reservoirs for flood control are presented in Table 2 and tested in this paper. Some control decisions are made for these types of heuristic controls, such as the critical water surface depth to open the valves for the on/off approach and the required retention time after a storm event to start releasing water for the detention control approach. These parameters were assumed as hcr=3  m (i.e., representing a critical level in terms of flood management) and td=6  h (i.e., representing a minimum detention time for sedimentation).Table 2. Types of static rule-based controlsTable 2. Types of static rule-based controlsTypeDescriptionPassive controlValve always 100% openDetention controlIf an event occurs, valve opening = 0%After the event, valve opening = 0% for tdElse, valve opening = 100%On/offIf hhclim))Δtwhere hlimc = maximum allowable water surface depth in the channel; ℭ = set of assessed controls; m = assessed control; and T = duration set of the analysis. Other similar metrics applied to real-time controls can be found in Wong and Kerkez (2018), Shishegar et al. (2019), and Schmitt et al. (2020).Results and DiscussionThis section presents results and discussions on how MPC can increase flood performance compared to other reactive control strategies. First, we assess the control performance to critical consecutive design storms in San Antonio in the section “Design Storms.” Following that, we perform a continuous simulation with observed data and assess the performance of MPC in contrast to reactive controls in the section “Continuous Simulation.”Design StormsThe comparison between reactive and predictive control approaches is shown in Fig. 6. The events produced similar peaks because of the initial saturation provided from the first storm. Outflow peaks in the reservoir increased from the first to the second storm, even though the runoff from the catchment decreased. The performance summary of rule-based control against model predictive control is given in Table 4. The control approach that generated the highest outflow peak was the detention control. This approach basically started to release the stored water at 18 h from the beginning of the first event (i.e., 6 h from the end of the first storm event) and it matched with the start of the intense inflow runoff volume at the detention basin originated from the second storm. Only the detention control and the on/off strategy reached the spillway elevation and hence had the higher outflows. For the MPC control, smaller flows than all RBCs were observed. From these results it is noted that reactive controls, especially the ones primarily designed for water quality purposes, can increase flood and consequently erosion downstream for large storm events (Schmitt et al. 2020).Table 4. Comparison within control strategies for a consecutive 25-year, 12-h, and 10-year, 12-h design storm in San Antonio for a watershed in a recharge zoneTable 4. Comparison within control strategies for a consecutive 25-year, 12-h, and 10-year, 12-h design storm in San Antonio for a watershed in a recharge zoneType of controlControl strategyPeak flow reduction first peak (%)Peak flow reduction second peak (%)Relative maximum flood depth (%)Relative control effort (%)Flood duration (h)Static/reactivePassive73.0868.6726.570.010.40On/off37.9727.81148.8993.4610.98Detention control15.011.09229.8715.588.38DLQI72.1167.5430.1861.3711.15DLQR72.4267.9129.0153.7410.83Optimization-based/predictiveInterior-point optimizer79.3976.022.5548.912.2The detention pond dumps faster on the passive control strategy as expected. However, the temporal evolution of the water surface depth in the reservoir had a similar behavior compared to DLQI and DLQR controls. This suggests that using reactive controls for relatively large floods might be as effective as the passive control. For the DLQI it is noted that the state hr(k) is steered to the reference setpoint (r(k)=1  m) where the system dynamics are not changing dramatically (i.e., cases where the inflow is steady). This type of control strategy is more suitable where constraints of minimum water depths are required in the reservoir (e.g., wetlands). Significant differences between the control approaches, however, are depicted in the channel stage variation over time. For all reactive controls, periods of flooding occurred, whereas the MPC control largely avoided it, although we assumed a perfect 8-h rainfall forecasting.For the MPC, a solution of a nonlinear optimization problem was developed for each control horizon. We used the interior-point method, which is a gradient-based method dependent on the initial estimated control decision. Therefore, to initialize the optimization algorithm, we assume this initial point as a random combination from the previous control signal. Basically, a normal distributed random number with 0 average and 0.2 variance is generated (δ=∼N(0,0.2)) and applied to the previous control, such that the initial point for the optimization is u0(k+1)=u(k)(1+δ). Defining this initial estimate, however, does not imply solutions only within this interval but actually can substantially decrease the computational time by starting in initial points near the previous control state.Continuous SimulationThe comparison between reactive approaches with an optimized-based/predictive control approach in a continuous simulation is shown in Fig. 7 and the trade-offs between control effort and flood duration for the continuous modeling are shown in Fig. 8. In this analysis, we assume a 2-h rainfall forecasting, showing the performance of RTC for low degrees of uncertainty in rainfall. In Fig. 7, a detail of the storm event that occurred in May 1, 2021, shows the hydrographs and the water surface depths in the reservoir and channels as well as the control schedules over time.This rainfall event was rapid and intense, producing a high inflow peak because of the initial saturation of the soil from previous rainfall events. The hydrograph shows that the on/off and detention control had the highest outflows peak. The valves were mostly closed for all reactive control strategies in the arriving of the inflow peak because no relatively high-water stages were yet observed. However, the on/off started this event with approximately 3 m of water stage (i.e., its control strategy is regulated by this water depth).While almost all other control strategies (i.e., except detention control) were able to release the accumulated volume with reasonable maximum channel water depths, the on/off strategy produced the highest outflow peak. It occurred because of opening both orifices late without avoiding spillway outflow. The risk of flooding is generally increased with this strategy. Although the detention control had some of the inflow volume spilled, all the stored volume corresponding to approximately 2.5 m (i.e., 5.5–3.0  m) of the pond stored volume [Fig. 7(b)] was successfully released 6 h after the end of the rainfall event. As shown in the section “Design Storms,” this available time might not always be feasible if a new storm event occurs in this interval.The DLQI, DLQR, and passive control strategies had similar outflows. The common reaction of the regulators is to stabilize the changing dynamics described by ψ(k,x(k),u(k)) due to qoutw(k) (Fig. 3); therefore, rapid valve openings is a common solution chosen by these controls depending on the tuned matrices Q and R. One also can tune these matrices differently to favorably change the relative importance of variations in control schedule and water stage over time.It is noted from Fig. 7 that the the detention and on/off controls had the highest stored water surface depths, which is in accordance with Sharior et al. (2019). The latter control, however, returned for the regulated water surface depth relatively quickly because of the decision to open the valves after reaching 3 m, whereas the detention control only released the flow hours afterward. Moreover, even after the event, some volume is still stored for all controls. Outflow discharge only occurs when hr(k)>h0+hm×Dh, where Dh is the hydraulic diameter. Thus, the minimum water surface depth is approximately 0.24 m.An interesting result is the fact that the predictive control had a relatively larger stored water surface depth in the reservoir compared to the other controls. Despite this fact, it still regulated the maximum water surface depth in the channel below the reference href of 1.8 m, as shown in Fig. 7(c). During the peak of the storms, the predictive control decided to partially close the valves, which was a feasible decision because the water surface depth in the pond was below the spillway elevation and no intense inflow was expected within the next few hours. This type of decision is counterintuitive and pinpoints the benefit of using predictive optimization-based control algorithms rather than solely expertise-based manual operations in real time.Frequency plots in Fig. 9 show the exceedance probability of outflows and stages in the reservoir and channel. The on/off control is more likely to release high flows, which can cause erosion (Schmitt et al. 2020). This control is similar to a pond with the spillway in an elevation of 2.5 m above the ground. It certainly has the benefits of water quality enhancement because of larger retention times (Sharior et al. 2019; Wong and Kerkez 2018), but this approach was the one with higher outflows and hence water surface depths in the channel. The DLQR had the more spread-out flow–duration curve, indicating that this control typically can release flows throughout a larger period of time in a relatively small magnitude. This is particularly important because the pond can be slowly emptied when no future inflow is forecasted. The DLQI, in addition to the on/off control, had a flat duration curve around their regulated water surface depths of 3 and 1 m (Table 2), respectively, which indicates that they can satisfy their control goals over a rainy season.Results in Figs. 9(c and d) show an estimate of the flood probability in the channel. They indicate that the passive control is the most suitable reactive control to avoid flood in the channel. However, the DLQR had approximately the same effect for mitigating the flood in the channel, but added the benefit of maintaining water in the channel and in the reservoir for longer periods of time. This type of RTC can therefore increase the period of time in the channel with some sort of base flow (Xu et al. 2021). The predictive control had no probability of flood in the channel, as desired from the optimization function.Another interesting analysis is the trade-offs between the control efforts and the flood duration. The largest control efforts were provided by the on/off control, detention control, and DLQI, as shown in Fig. 8. Because we tune Q with a high weight in the output deviation (e˙=r(k)−y(k)), the control law rapidly tries to steer the system back to the defined reference setpoint, thus generating intense control efforts. This tuning can represent the operator preferences (Troutman et al. 2020; Fraternali et al. 2012) and could be a time-varying function in terms of the rainfall forecastings, although most control algorithms assume a constant value. The control approach with least control effort was the predictive control, which nearly had 10% of the detention control approach. The MPC problem for the continuous simulation was solved in approximately 30 min in a single core machine.Conclusions, Limitations, and Future WorksThe real-time control of watersheds, reservoirs, and channels can decrease the risk of flooding in stormwater systems by better controlling their flow release over time. The application of different types of optimization-based and reactive controls algorithms was tested in a case study, and given the questions posed in the section “Mathematical Model Application,” the results support the following conclusions: 1.The reactive controls have lower flood mitigation performance when compared to the predictive strategy, especially the water quality controls of the on/off control and the detention control.2.The MPC control outperforms all other control strategies for flood mitigation performance and also requires less control effort.3.The DLQR and DLQI are as effective as the passive control for flood mitigation in the channel. However, they add more specific benefits such as maintaining a desired water surface depth in the pond (DLQI) and increasing residence times and low flows (DLQR). However, tuning the weighing matrices Q and R from the DLQR and LQR optimization problems can be complex because of different units in the objective function. Therefore, using normalized objective functions could be an alternative.4.Although the reactive controls were each applied for a control interval of 15 min and the predictive control for a control interval of 60 min, the larger opportunity to change the valve operation did not compensate the lack of predictability of the states for the reactive controls for either large critical design storms events or for observed rainfall events.Although only flood performance is evaluated in this paper, this type of modeling approach can be flexible enough to combine different control purposes according to the expected rainfall forecastings. Depending on the estimated future states of the system, one can change the objective functions of the system to shift from a flood-focused strategy to a water quality control, increasing detention times. Moreover, because uncertainty in forecasting increases over the predicted horizon, it is possible to tune the weight matrices in the objective function following the uncertainty associated in the forecasting and thus giving more importance for short-duration forecastings by increasing the cost weights associated with this duration.Future application of the developed approach will consider scenarios of uncertainty in rainfall, model parameters, and measurements noise. Moreover, the approach needs to be tested in real case studies, such as the Upper San Antonio River watershed, which contains some of the most advanced flood protection systems in the United States. Comparisons of the simplified dynamical system with the state-of-the-art hydraulic and hydrologic models, such as SWMM, GSSHA, HEC-RAS, InfoWorks, and many others, are also warranted. This methodology can be easily expanded for systems with many watersheds, reservoirs, and channels.Data Availability StatementSome or all data, models, or code generated or used during the study are available in a repository or online in accordance with funder data retention policies. All functions, scripts, and input data including the state-space nonlinear model, the model predictive control algorithm, and the rule-based algorithms are available in Gomes (2021).AcknowledgmentsThe authors gratefully acknowledge the support by the City of San Antonio, the San Antonio River Authority, and the National Science Foundation (NSF) under Grant 2015671 and CAPES/PROEX-Coord. Aperf. Pessoal Ensino Superior, Programa de Excelência at Graduate Program of Hydraulics and Sanitation (PPGSHS) at Escola de Engenharia de São Carlos, University of São Paulo, Brazil.References Akan, O. A. 1993. Urban stormwater hydrology: A guide to engineering calculations. Boca Raton, FL: CRC Press. Brasil, J., M. Macedo, C. Lago, T. Oliveira, M. Júnior, T. Oliveira, and E. Mendiondo. 2021. “Nature-based solutions and real-time control: Challenges and opportunities.” Water 13 (5): 651. Duchesne, S., A. Mailhot, E. Dequidt, and J. P. Villeneuve. 2001. “Mathematical modeling of sewers under surcharge for real time control of combined sewer overflows.” Urban Water 3 (4): 241–252. Fraternali, P., A. Castelletti, R. Soncini-Sessa, C. V. Ruiz, and A. E. Rizzoli. 2012. “Putting humans in the loop: Social computing for water resources management.” Environ. Modell. Software 37 (Nov): 68–77. Fry, T. J., and R. M. Maxwell. 2018. “Using a distributed hydrologic model to improve the green infrastructure parameterization used in a lumped model.” Water 10 (12): 1756. Furl, C., D. Ghebreyesus, and H. O. Sharif. 2018. “Assessment of the performance of satellite-based precipitation products for flood events across diverse spatial scales using GSSHA modeling system.” Geosciences 8 (6): 191. García, L., J. Barreiro-Gomez, E. Escobar, D. Téllez, N. Quijano, and C. Ocampo-Martinez. 2015. “Modeling and real-time control of urban drainage systems: A review.” Adv. Water Resour. 85 (Nov): 120–132. Gasper, R., A. Blohm, and M. Ruth. 2011. “Social and economic impacts of climate change on the urban environment.” Curr. Opinion Environ. Sustainability 3 (3): 150–157. Giacomoni, M. H., and J. Joseph. 2017. “Multi-objective evolutionary optimization and Monte Carlo simulation for placement of low impact development in the catchment scale.” J. Water Resour. Plann. Manage. 143 (9): 04017053. Gomes, M. N., Jr., M. H. Giacomoni, and E. M. Mendiondo. 2021a. “Simpósio Brasileiro de Recursos.” In The role of raster resolution into overland flow and total suspended solids modeling in small urban catchments. Belo Horizonte, Brazil: Associação Brasileira de Recursos Hídricos. Gomes, M. N., Jr., M. H. Giacomoni, A. T. Papagiannakis, E. M. Mendiondo, and F. Dornelles. 2021b. “Spatial assessment of overland flow, pollutant concentration, and first flush using a 2-D non-point source pollution and hydrological model for urban catchments.” In Proc., World Environmental and Water Resources Congress 2021: Planning a Resilient Future along America’s Freshwaters—Selected Papers from the World Environmental and Water Resources Congress 2021, 397–413. Reston, VA: ASCE. Horváth, K., M. Petreczky, L. Rajaoarisoa, E. Duviella, and K. Chuquet. 2014. “MPC control of water level in a navigation canal—The Cuinchy-Fontinettes case study.” In Proc., 2014 European Control Conference (ECC), 1337–1342. New York: IEEE. Joseph-Duran, B., C. Ocampo-Martinez, and G. Cembrano. 2015. “Output-feedback control of combined sewer networks through receding horizon control with moving horizon estimation.” Water Resour. Res. 51 (10): 8129–8145. Kearfott, R. B. 2013. Vol. 13 of Rigorous global search: Continuous problems. New York: Springer. Kessler, R. 2011. “Stormwater strategies: Cities prepare aging infrastructure for climate change.” Environ. Health Perspect. 119 (12): A514–A519. Kollet, S. J., and R. M. Maxwell. 2006. “Integrated surface-groundwater flow modeling: A free-surface overland flow boundary condition in a parallel groundwater flow model.” Adv. Water Resour. 29 (7): 945–958. Li, J. 2020. “A data-driven improved fuzzy logic control optimization-simulation tool for reducing flooding volume at downstream urban drainage systems.” Sci. Total Environ. 732 (Aug): 138931. Lund, N. S. V., A. K. V. Falk, M. Borup, H. Madsen, and P. Steen Mikkelsen. 2018. “Model predictive control of urban drainage systems: A review and perspective towards smart real-time water management.” Cri. Rev. Environ. Sci. Technol. 48 (3): 279–339. Maxwell, R. M., et al. 2014. “Surface-subsurface model intercomparison: A first set of benchmark results to diagnose integrated hydrology and feedbacks.” Water Resour. Res. 50 (2): 1531–1549. McCarthy, G. T. 1938. “The unit hydrograph and flood routing.” In Proc., Conf. of North Atlantic Division, 608–609. New London, CT: USACE. NOAA (National Oceanic and Atmospheric Administration). 2021 Climate data online. Washington, DC: NOAA. Ocampo-Martínez, C., and V. Puig. 2010. “Piece-wise linear functions-based model predictive control of large-scale sewage systems.” IET Control Theory Appl. 4 (9): 1581–1593. Pappas, T., A. Laub, and N. Sandell. 1980. “On the numerical solution of the discrete-time algebraic Riccati equation.” IEEE Trans. Autom. Control 25 (4): 631–641. Pappenberger, F., K. Beven, M. Horritt, and S. Blazkova. 2002. “Uncertainty in the calibration of effective roughness parameters in HEC-RAS using inundation and downstream level observations.” J. Hydrol. 302 (1–4): 46–69. Rossman, L. A. 2010. Storm water management model user’s manual, version 5.0. Cincinnati: National Risk Management Research Laboratory, Office Of Research And Development US Environmental Protection Agency. Schmitt, Z. K., C. C. Hodges, and R. L. Dymond. 2020. “Simulation and assessment of long-term stormwater basin performance under real-time control retrofits.” Urban Water J. 17 (5): 467–480. Sharior, S., W. McDonald, and A. J. Parolari. 2019. “Improved reliability of stormwater detention basin performance through water quality data-informed real-time control.” J. Hydrol. 573 (Jun): 422–431. Shishegar, S., S. Duchesne, and G. Pelletier. 2019. “An integrated optimization and rule-based approach for predictive real time control of urban stormwater management systems.” J. Hydrol. 577 (Oct): 124000. Te Chow, V. 2010. Applied hydrology. New York: Tata McGraw-Hill Education. Troutman, S. C., N. G. Love, and B. Kerkez. 2020. “Balancing water quality and flows in combined sewer systems using real-time control.” Environ. Sci.: Water Res. Technol. 6 (5): 1357–1369. Vose, M. D. 1999. The simple genetic algorithm: Foundations and theory. Cambridge, MA: MIT Press. Wang, S., A. F. Taha, and A. A. Abokifa. 2021. “How effective is model predictive control in real-time water quality regulation? State-space modeling and scalable control.” Water Resour. Res. 57 (5): e2020WR027771. Wing, O. E., N. Pinter, P. D. Bates, and C. Kousky. 2020. “New insights into us flood vulnerability revealed from flood insurance big data.” Nat. Commun. 11 (1): 1–10. Wong, B. P., and B. Kerkez. 2018. “Real-time control of urban headwater catchments through linear feedback: Performance, analysis, and site selection.” Water Resour. Res. 54 (10): 7309–7330. Xu, W. D., M. J. Burns, F. Cherqui, and T. D. Fletcher. 2021. “Enhancing stormwater control measures using real-time control technology: A review.” Urban Water J. 18 (2): 101–114. Zhang, W., G. Villarini, G. A. Vecchi, and J. A. Smith. 2018. “Urbanization exacerbated the rainfall and flooding caused by Hurricane Harvey in Houston.” Nature 563 (7731): 384–388.

Source link

Leave a Reply

Your email address will not be published. Required fields are marked *