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 h
Source link
