Decentralized Control for Urban Drainage Systems using Replicator Dynamics

This paper proposes a decentralized control scheme that mitigates floods in urban drainage systems (UDSs). First, we develop a partitioning algorithm of the UDS relying on a graph model of the system. Once this is done, we design a local controller for each partition based on the replicator dynamics model (a set of differential equations that describes the evolution of a population of players involved in a strategic game). The decentralized nature of the proposed strategy makes it suitable for applying it in large-scale systems. Stability of the closed-loop system is proved by using Lyapunov theory. Furthermore, we simulate the performance of the decentralized control scheme in two case studies. One of them models part of the Bogotá (Colombia) stormwater UDS. Finally, we compare the proposed technique with two widely used methods for-real time control of UDSs, i.e., constrained linear quadratic regulator (LQR) and model predictive control (MPC).


I. INTRODUCTION
Urban drainage systems (UDSs) collect, store, and transport residual water to wastewater treatment plants and water bodies. UDSs are composed of an arrangement of channels which are connected by inspection and collection chambers. The wastewater flow along UDS is regulated by using measuring equipment (e.g., level and flow sensors) and actuators (e.g., valves and gates).
Normally, UDSs are designed for critical rain and demand scenarios up to certain magnitude [1]. However, in some situations, the wastewater flow that is conveyed by the UDS surpasses its capacity. For this reason, the system becomes saturated and flooding happens. In fact, some recent studies have shown that flood risk is increasing, specially in large cities [2], due to several factors, e.g., an increment of magnitude of extreme precipitation as a result of global warming, a rapid growth of population that stresses drainage infrastructures, and construction of civil infrastructure that has reduced the soil capacity to absorb rainwater [3]. Floods have severe impacts including loss of human life, destruction of property, and health problems.
For the reasons described above, flooding mitigation has gained broad research attention. There are two main approaches: the first one consists in upgrading UDSs infrastructure [4] (e.g., by the construction of new sewer collectors or re-sizing of the existing ones). The second approach is to use real-time control techniques that take advantage of the system's actuators to modify flow conditions throughout the UDS. Implementation of real-time controllers is appealing because it is more cost effective than enlarging existing drainage infrastructures. Therefore, several control techniques have been proposed [5]- [9]. One common approach is to develop simulation tools of UDSs and employing them as black-box models in heuristic-based optimization algorithms [10]- [13]. However, these algorithms lack a rigorous analysis that guarantees, for example, their stability. Another strategy relies on optimal control. For instance, [14] employs a multivariable linear quadratic regulator (LQR) that uses all the available storage capacity of the system and maximizes its outflow to treatment plants. Other authors have also used model predictive control (MPC) for UDSs management [15]- [19].
A widespread problem of most of the aforementioned techniques is that they require the design of a central authority in charge of the whole UDS. This centralized controller is responsible for handling all the actuators and processing information from all system's sensors. Since the size of UDSs is generally of large-scale nature, implementing a centralized controller is problematic. Instead, it is desirable to design decentralized control schemes, i.e., schemes with multiple controllers, where each one of them has partial information of the UDS. Decentralized control of UDS has also been tackled in the literature. For instance, using rule based controllers, as in [20], where fuzzy logic is employed to manage the wastewater level of storage units. Nevertheless, this approach does not consider any coordination of storage units, which largely reduces the benefits of the strategy [21]. On the other hand, some authors have used distributed versions of model predictive control to overcome the requirement of a centralized authority (e.g., see [22]). However, in these schemes the computation of control signals requires high-performance processors because of the complexity of existing algorithms.
Recently, game theory and population dynamics based controllers have gained relevance due to their low computational burden and simplicity for designing decoupled schemes using them (e.g., see [23]- [29]). In fact, there are multiple engineering problems that have been addressed by means of game theory and population dynamics. For instance, coordination of robot networks [30], wind farm optimization [31], demand response in electrical grids [32], traffic assignment [33], charging of large populations of electric vehicles [34], control of epidemics [35], and so forth. It is important to point out that population dynamics based controllers have been applied to water allocation problems. For instance, [36] poses a game between the water resource manufacturers and regulators to find an efficient water allocation in a river basin. In particular, population dynamics have also been used to control UDSs. As described in [37], this kind of control strategy has shown suitable performance for mitigating floods.
Population dynamics play the key role in the control strategy that we propose in this article. Specifically, we develop a decentralized architecture that divides the UDS in subsystems handled by local controllers. Each one of these controllers implements a population dynamics algorithm, which is designed to make an efficient wastewater allocation among the storing units of the sub-system. Efficiency of the proposed strategy is proved using Lyapunov analysis. Furthermore, we test the performance of the controlled UDS by means of simulations under two scenarios. One of them is a case study that uses the model of part of the Bogotá (Colombia) stormwater UDS.
It is worth mentioning that, compared to [23] and [37], where a similar control strategy is used, our paper has three main advantages: First, the formalization of a partitioning algorithm that covers a wider range of UDSs topologies. Second, the wastewater allocation process performed by each local controller takes into account the available capacity of all storing units of the handled sub-system (this is not the case of the strategy proposed in [23] and [37] where only upstream storing units are involved in the resource allocation process). Finally, we have removed avoidable constraints on the control signals that are often imposed when implementing population dynamics based algorithms (specifically, in population dynamics based controllers, the sum of control actions are generally required to remain constant along the time). These constraints unnecessarily reduce the degrees of freedom of control actions.
Regarding partitioning methods that divide UDSs into subsystems, some procedures have been reported in the literature. For instance, [38] proposes a cascade configuration of MPC controllers. First, the UDS is divided into subsystems from downstream to upstream reservoirs. Then, the control strategy starts by computing the outputs of controllers in charge of downstream sub-systems. Once this is done, information is sent to upstream sub-systems so that their controllers can compute the control signals. A drawback of this strategy is that parallelization is not possible since upstream controllers require the information of downstream controllers. Another partitioning method is the one described in [39], where the authors develop an algorithm based on the Strahler number. The idea is to decompose the UDS in hierarchically interconnected sub-systems to increase the system's resilience. This strategy is intended to be implemented at the planning stage of the UDS rather than in its operation stage. Hence, the resulting partitions do not fit well with realtime controllers. On the other hand, [40] develops a gossipbased algorithm that coordinates local controllers in a UDS divided into sub-systems. However, the partitioning method employed by the authors can only be applied to UDSs with tree topologies. Furthermore, [41] devises graph partitioning algorithms to group sections of a UDS into semi-distributed sub-systems. Although these algorithms perform well to accelerate accurate simulations of UDSs, they are not designed for control purposes. Other authors use a geographical proximity criterion for UDSs partitioning [42], which minimizes the cost of the communication infrastructure required by each local controller. Nevertheless, this criterion does not take into account the information requirements of the controllers. Hence, the performance of the closed-loop system may worsen. Another partitioning paradigm that emerges in optimal control strategies for UDSs consists in dividing the optimization problem into sub-problems (which are related to sub-systems of the UDS) [22]. Each sub-problem is solved by a local controller while all the solutions are coordinated by means of a distributed algorithm (e.g., the alternating direction method of multipliers (ADMM) [43]). The main drawback of this paradigm is that coordination algorithms require a communication network that connects all subsystems (generally, such networks are modeled as connected graphs). In addition, convergence rates of distributed methods are quite lower compared to centralized approaches. Finally, some papers also propose to separate the collection of wastewater and surface run-off at the planning stage (see e.g., [44]). This division alleviates the risk of overflows. However, separate drainage systems require higher infrastructure and maintenance costs [45]. Beyond the drawbacks that we have pointed out for each of the cited partitioning methods, the key problem is that none of these methods produces subsystems with topologies that fit the controller proposed in this article. Thus, our partitioning algorithm is a central piece in the design of the management strategy presented in this paper. We want to highlight that our partitioning algorithm is oriented to decompose the UDS into sub-systems that use a local control policy, where each sub-system is capable to compute the control signals independently of the state of the other sub-systems. This fact offers inherent resilience to failures (which is typical of decentralized control schemes) because the malfunction of a subsystem does not compromise the performance of the entire UDS.
In addition, the fact that the control strategy developed in this paper uses population dynamics leads to advantages compared to other decentralized methods. On the one hand, population dynamics only requires the evaluation of a set of algebraic expressions to update the control actions. Therefore, compared to optimization-based techniques, our approach avoids numerical issues raised by the employment of optimization solvers. Indeed, in Section VII, we show by means of simulations that the strategy proposed in this article outperforms some optimization-based approaches and decreases the computational burden because of the different control paradigm employed. On the other hand, the analytical nature of population dynamics allows us to perform mathematical analysis for guaranteeing key features of the closed-loop controlled system such as stability. This feature is missing in some UDS control methods, especially in those based on heuristics.
The remainder of this paper is organized as follows: Section II covers notation and some general concepts about population dynamics. In Section III, we describe the problem statement and the model of the UDS employed for developing the control strategy. Sections IV and V show the partitioning algorithm and the design of local controllers using population dynamics, respectively. In addition, a stability analysis of the controlled system is presented in Section VI. Furthermore, Section VII discusses the simulation results obtained when the proposed strategy is implemented in two case studies. This section also presents comparisons with commonly used methods for real-time control of UDSs. Specifically, we compare our method with two optimizationbased controllers, i.e., constrained LQR and MPC. Finally, some conclusions are drawn in Section VIII.

A. NOTATION
Throughout this article, bold and calligraphic letters denote vectors and sets, respectively. Variables that depend on time are explicitly written as functions of t, e.g., x(t). As usual, the cartesian product of two sets A and B is denoted by A × B, and the number of elements of the set A is denoted by |A|.
Nonnegative and positive real numbers are denoted by R ≥0 and R >0 , respectively. Finally, R n ≥0 is the set of vectors of size n whose entries belong to R ≥0 .

B. REPLICATOR DYNAMICS MODEL
The problem of flooding can be mitigated if wastewater is efficiently distributed in the storing units (chambers) of the UDSs. Consequently, an appropriate wastewater allocation mechanism is required. If we consider the sewage in the UDS as a resource, traditional resource allocation methods (such as those based on static optimization) are not suitable to address the problem under consideration because the resource is not static. Indeed, UDS is continuously receiving (e.g., from rain) and discharging sewage. Thus, a dynamic resource allocation mechanism, as the case of the replicator dynamics, is required. Broadly speaking, replicator dynamics describe how a mass of individuals (e.g., animals, agents) tries to efficiently distribute itself among different strategies (e.g., habitats, actions). Efficiency of the distribution is measured by means of payoff functions known as fitness functions.
For a better understanding of the role that replicator dynamics play in wastewater allocation, let us consider the following analogy. Assume that the mass of individuals is the wastewater of the UDS and the strategies are the chambers of the system. Under these assumptions, the key idea is to design proper fitness functions in such a way that the replicator dynamics allocate the sewage (mass of individuals) in emptier chambers (more efficient strategies).
Since the control strategy proposed in this article is based on the replicator dynamics, it is necessary to describe its mathematical model (introduced by Taylor and Jonker in 1978 [46]). The replicator equation is a classic model that captures the natural selection process in a population of individuals that can choose between different habitats or strategies. The key elements of the model are: A set of available strategies, which is denoted by S = {1, . . . , s} (in this case we have s strategies); a population mass P ; a set of variables p 1 (t), . . . , p s (t), where p i (t) describes the portion of individuals choosing the i-th strategy at time t; and a population state p(t) = [p 1 (t), . . . , p s (t)] ⊤ . Formally, the replicator dynamics model is given by the following differential equation: where f i (t) is a real-valued fitness function, which captures the payoff perceived by the individuals playing the i-th strategy,f (t) = 1 P j∈S p j (t)f j (t) is a weighted average fitness, and β ∈ R >0 is the population growth rate. Notice that according to (1), populations of most successful individuals (i.e., individuals that play strategies with higher-than-theaverage fitness functions) tend to grow, while least successful populations decrease. This behavior is closely related with the biological mechanism of natural selection, i.e., survivalof-the-fittest. One important property of this model is next described in Lemma 1. The importance of Lemma 1 lies in the fact that it guarantees the evolution of the population state p(t) in the nonnegative orthant provided that p(0) ∈ R n ≥0 . Furthermore, the population mass P is preserved along the time. These two properties of p(t) will be exploited in the controller design.
Another key feature of the replicator equation is that the fitness functions of non-extinct strategies 1 reach the same value at equilibrium. Formally, if p * i , p * j ∈ R >0 , then f * i = f * j , for all i, j ∈ S, where p * i and f * i denote the values of p i (t) and f i (t) in steady-state, respectively. This fact motivates the use of replicator dynamics in problems where the objective is to attain consensus with respect to a certain variable. For instance, the temperature of the rooms in a building [48], the frequency of electric generators [49], the illuminance level in offices with several lamps [23], or the remaining storage capacity in multi-tank systems, which is the case of this article.

Including Populations Upper Bounds
The population of individuals playing each one of the available strategies is lower-bounded by zero under the replicator dynamics in (1) (cf., Lemma 1). Nevertheless, this portion is not bounded above by any quantity lower than the population mass P . Although upper bounds are not considered in (1), they can be included by adding in the formulation of the replicator equation a term related to the maximum size of each population, i.e., a term that prevents p i (t) to exceed an upper boundp i . The described modification, which has been explored in [50], is specifically done in the fitness functions as follows: wheref i (t) is the payoff associated with the i-th strategy without considering upper bound, ε ∈ R ≥0 is a small nonnegative scale factor, and b i (p i (t),p i ) is a function that prevents p i (t) to exceed the upper boundp i . Modified fitness function has the following characteristics: R is a Lipschitz continuous and monotonically decreasing function with respect to its first argument.
The behavior of the population playing the i-th strategy (i.e., p i (t)) under the fitness functions in (2) can be analyzed considering two cases. The first case is when the variable p i (t) is far fromp i . Notice that it is possible to choose a sufficient small value for ε such that the effect of b i (p i (t),p i ) on the fitness f i (t) is negligible if p i (t) is far fromp i . In this case, the fitness function is quite close to the payoff associated with the i-th strategy without considering an upper bound, i.e., f i (t) ≈f i (t) according to (2). The second case is when p i (t) tends to the valuep i . Under this scenario, although the value of ε is small, the function εb i (p i (t),p i ) takes an arbitrarily large negative value due to P2. Therefore, the fitness f i (t) → −∞. Notice that since f i (t) → −∞, the payoff perceived by the i-th population is lower than the average, i.e., f i (t) <f (t). Hence,ṗ i (t) < 0 cf., (1) , and thus p i (t) does not exceed its corresponding upper bound p i . Summarizing, p i (t) ∈ [0,p i ). The previous discussion is illustrated in Figure 1 value of the barrier function

A. PROBLEM STATEMENT AND OVERVIEW OF THE PROPOSED SOLUTION
We address the problem of overflooding in UDSs. This problem occurs when the wastewater exceeds the storage capacity of any of the UDS's collectors. Therefore, overflooding can be avoided if wastewater is efficiently allocated among all the UDS's collectors by handling the system's actuators (valves in our case). Reaching an efficient wastewater allocation is the key idea of the proposed control strategy. Specifically, we seek that all collectors have the same remaining capacity for storing wastewater. In other words, our objective is trying that all UDS's collectors have the same empty space during mostly of the time. If this happens, no collector will be overloaded compared to the other ones. Therefore, no flooding occurs provided that there is available storage capacity in the UDS. On the other hand, UDSs are of large scale nature. Hence, decentralized control strategies suit best for these kinds of systems, i.e., we are looking to design a set of controllers, where each controller is in charge of part of the UDS. A constraint that must be satisfied by local controllers is that they are only able to use information from the subsystem they are handling. For attaining the control objective while meeting the local information constraint, we divide our strategy in two main stages: first, the UDS is partitioned in sub-systems that have some desired characteristics, and then a local controller based on population dynamics is designed for each sub-system obtained in the partitioning stage.
The partitioning stage uses a graph model of the UDS while the controllers design is based on a dynamic model derived from the mass conservation principle. Let us describe these two models.

B. GRAPH MODEL OF A UDS
Using the ideas in [23], we represent each collector of the UDS as a reservoir. A UDS with m reservoirs is modeled as a connected and directed graph G = (V, E), where V = {1, . . . , m} is the set of nodes of the graph, and E ⊆ V × V is the set of edges. In our case, the nodes are the reservoirs of UDS, and the edges represent connections between these reservoirs, e.g., the edge (i, j) ∈ E models that there is a connection allowing wastewater flow from the i-th reservoir to the j-th reservoir. We also assume that each edge has a valve that controls wastewater flow from one reservoir to another.
The UDS can be divided into sub-systems. Each subsystem belongs to one of the following two topologies.

Definition 1. (Divergent topology)
The divergent topology involves n + 1 reservoirs, and represents a situation in which there are n > 1 receptor reservoirs and only one source reservoir. Consequently, in a divergent topology there are n manipulated flows diverging from a unique source reservoir. Figure 2(a) shows a scheme representing such a topology. Formally, a sub-system with divergent topology can be represented as a graph where the nodes are the source reservoir and the receptor reservoirs, and the edges are the links between the source and receptor reservoirs. ♢

Definition 2. (Convergent topology)
The convergent topology involves n + 1 reservoirs, and represents a situation in which there are n ≥ 1 source reservoirs and only one receptor reservoir. Consequently, in a convergent topology there are n manipulated flows converging to a unique receptor reservoir. Figure 2(b) shows a convergent topology. Formally, a sub-system with divergent topology can be represented as a graph where the nodes are the receptor reservoir and the source reservoirs, and the edges are the links between the source and receptor reservoirs. ♢ Remark 1. If one sub-system has only one source reservoir and one receptor reservoir (see Figure 2(c)), then it can be represented as a convergent topology. Notice that divergent topologies require at least two receptor reservoirs.

C. DYNAMIC MODEL OF RESERVOIRS
We use the Muskingum model [51], which is based on the principle of mass conservation (this principle is widely used in the literature for control-oriented models [52]). Specifically, we represent each channel of the UDS as a reservoir containing a wastewater volume v i (t). Besides, we assume that each channel has a set of output gates whose opening percentage can be controlled. Dynamics of the wastewater volume of the i-th reservoir can be described by the following differential equation: where q in,i (t) and q out,i (t) are the reservoir inflow and outflow, respectively. Modeling of inflow and outflow depends on the topology of the sub-system to which the reservoir belongs. In Section V, we delve into the dynamic model of reservoirs depending on sub-systems topologies.
It is worth noting that convergent topologies have been addressed before in [23] and [37]. However, these papers do not consider the possibility that a UDS has sub-systems with divergent topology. Therefore, our approach covers a wider range of UDSs architectures.

IV. SYSTEM PARTITIONING
Partitioning consists in dividing the whole UDS in subsystems belonging to the two topologies described before. Besides, partitioning must satisfy the following constraints: C1. All resulting sub-systems must belong to one and only one topology (convergent or divergent). This is required because we only design two types of controllers, one per each topology. C2. All reservoirs of the UDS must belong to at least one sub-system. If one reservoir does not belong to any subsystem, its wastewater volume would not be controlled. C3. Each edge of the graph that models the UDS must belong to one and only one sub-system. This constraint prevents ambiguity in control actions as explained in Section V. Any partitioning algorithm satisfying these three constraints is suitable for the splitting of UDS. In particular, we propose an algorithm consisting in three main steps: 1) identification of divergent topologies; 2) removal of the edges that belong to divergent topologies; and 3) identification of convergent topologies.
Before stating the partitioning algorithm, let us define two concepts related to directed graphs. Definition 3. In-neighborhood of a node: Let G = (V, E) be a directed graph. The in-neighborhood of the node i ∈ V, which is denoted by N in i (G), is the set of nodes satisfying the following property: Notice that if G is a graph that models a UDS, then N in i (G) is the collection of all reservoirs having an outflow that goes to the i-th reservoir. Definition 4. Out-neighborhood of a node: Let G = (V, E) be a directed graph. The out-neighborhood of the node i ∈ V, which is denoted by N out i (G), is the set of nodes satisfying the following property: Notice that N out i (G) is the collection of all reservoirs having an inflow that comes from the i-th reservoir.
Using these definitions, we describe the three stages of our partitioning algorithm.

A. STEP 1 (IDENTIFICATION OF DIVERGENT TOPOLOGIES)
According to Definition 1, a sub-system having divergent topology has one source node and several (at least two) receptor nodes. Hence, for identifying all the UDS subsystems with divergent topology, it is enough to find all the source nodes having more than one out-neighbor.
Formally, the set of source nodes of divergent topologies, which we denote by V s , is The receptor nodes associated with node i ∈ V s are N out i (G). Thus, a sub-system with divergent topology associated to the source node i ∈ V s is the graph G div are the source node and its receptors, i.e., Besides, the edges of G div i are the connections between the source node and its receptors, i.e.,

B. STEP 2 (REMOVAL OF THE EDGES BELONGING TO DIVERGENT TOPOLOGIES)
From the graph that models the UDS, we create a new graph subtracting all edges belonging to divergent topologies. Formally, the new graph is G ′ = (V, E ′ ), where the set of nodes is the same as the one of G, and the set of edges is

C. STEP 3 (IDENTIFICATION OF CONVERGENT TOPOLOGIES)
Similarly as in Step 1, we identify the sub-systems having convergent topology. However, we now use the new graph G ′ . We recall that a sub-system having divergent topology has one receptor node and at least one source node (see Definition 2). Hence, for identifying all the UDS subsystems with convergent topology, it is enough to find all the receptor nodes having at least one in-neighbor.
Formally, the set of receptor nodes of convergent topologies, which we denote by V r , is Hence, a sub-system with convergent topology associated to the receptor node i ∈ V r is the graph G cnv ). The nodes of G cnv i are the receptor node and its sources, i.e., V cnv Besides, the edges of G cnv i are the connections between the receptor node and its sources, i.e., E cnv i = {i} × N in i (G ′ ). All the partitioning algorithm steps described above are summarized in Figure 3.

D. PARTITIONING ALGORITHM EXAMPLE
To see how the partitioning algorithm works, we apply it in the UDS example shown in Figure 4(a), whose associated graph is depicted in Figure 4(b). The three steps of the partitioning algorithm are illustrated in Figure 5. Next, we describe each step. Figure 5(a) shows the result of applying the first step of the algorithm. All nodes having more than one out-neighbor (in this case, nodes 1 and 4) are source nodes of a divergent topology, i.e., V s = {1, 4}. Receptors of node 1 are nodes 3 and 4, and receptors of node 4 are nodes 5 and 6. Therefore, there are two sub-systems with divergent topology: G div 1 shown in blue, and G div 4 shown in red. Figure 5(b) shows the result of applying the second step of the algorithm. All edges of such sub-systems with diver-input G (graph that models the UDS) step 1.1 identifiy source nodes (nodes having more than one out-neighbor) identify receptor nodes (out-neighbors of each source node indentified in step 1.    shown in purple. Notice that, for this step, we use the graph G ′ , which is obtained after subtracting the edges belonging to divergent topologies from the graph that models the UDS. N out

V. DECENTRALIZED CONTROL OF UDS
Once the partition of the UDS is done, a controller is synthesized for each sub-system. As discussed in Section III-A, the control goal is to prevent overflows. Hence, wastewater must be efficiently allocated among the available reservoirs by controlling the system's valves. In our case, efficient wastewater allocation is achieved when the reservoirs of UDS have the same remaining capacity.

1) Dynamic model of sub-systems with flow divergent topology
As stated in Definition 1, sub-systems with flow divergent topology have one source reservoir and n receptor reservoirs (see Figure 2(a)). For convention, we use the indices 1, . . . , n, for identifying receptor reservoirs, while the source reservoir is marked with the index n + 1.
As explained in Section III-C, we represent each collector of the UDS as a reservoir of volumev i containing a wastewater volume v i (t). The model assumes that the outflows of a reservoir are proportional to its wastewater volume (for more details see [23]).
The differential equation describing the change along the time of the source reservoir's wastewater volume iṡ where k n+1 > 0 is a parameter that depends on the geometry of the source reservoir, x i (t) ∈ [0, 1] is the opening percentage of the valve that connects the source reservoir with the i-th receptor reservoir x i (t) = 0 means the valve is completely closed, and x i (t) = 1 means a completely opened valve , q in n+1 (t) ≥ 0 is the external inflow, and γ out n+1 (t) ≥ 0 is a coefficient related to the external outflow of the reservoir 2 .
The dynamic model of the receptor reservoirs' wastewater volume is given bẏ for all i = 1, . . . , n. Again, q in i (t) ≥ 0 is the external inflow, and γ out i (t) ≥ 0 is a coefficient related to the external outflow of the i-th receptor reservoir. Figure 6 shows a schematic representation of a sub-system with flow divergent topology. Notice that all the involved variables and parameters of the dynamic model described before are depicted in this figure.
v n+1 (t) q in n (t) Figure 6. A sub-system with divergent topology.
2) Replicator dynamics-based controller for flow divergent topology Based on the replicator dynamics model introduced in Section II-B, we design local controllers for each UDS sub-system having a flow divergent topology. Replicator dynamics-based controllers are fully characterized by a set of strategies, a population mass and a set of fitness functions. Let us define each one of these elements. In our game, each player has two options: increase the opening percentage of one of the tanks' valves, or decrease such a percentage for all the valves. If the player chooses the first option, then it must also select one of the n valves that connect the source tank with its n receptor tanks. Thus, we have n + 1 available strategies, i.e., the strategy set is S = {1, . . . , n + 1}.
Strategies 1, . . . , n, increase the opening percentage of the corresponding receptor tank, i.e., if a player selects the i-th strategy, where i = 1, . . . , n, then the opening percentage of the i-th receptor tank's valve increases. In this regard, the population playing this i-th strategy is exactly the opening percentage of the i-th receptor tank's valve, i.e., x i (t).
On the other hand, if a player adopts the strategy n + 1, it decreases the opening percentage of all valves. In this case, the population playing the strategy n + 1, i.e., x n+1 (t), does not have a physical meaning. However, it is equal to the sum of the remaining opening percentages of all valves.
According to previous definitions, the mass of the population involved in the replicator dynamics model is the sum of the opening capacity of the n valves that connect the source tank with its receptor tanks. Since we assume that x i (t) ∈ [0, 1], for all i = 1, . . . , n, then the population mass is equal to n.
Regarding fitness functions, the higher the remaining volume of a receptor tank, the more appealing to increase its wastewater inflow by increasing its valve opening percentage. This fact is captured in the following fitness functions: Notice that we have added a barrier term εb i (x i (t), 1) (as explained in Section II-B) to the first n fitness functions. In this barrier term, the upper bound of x i (t) is 1, because it is the maximum opening percentage. Notice also that the fitness function of the strategy n + 1 models the fact that if the remaining capacity of the source tank is high, then it is appealing to increase the wastewater stored in the source tank by decreasing the opening percentage of all valves.
Under previous assumptions, the replicator dynamicsbased controller is given by the following differential equations: where β is a tuning parameter of the controller, andf (t) = 1 n n+1 j=1 x j (t)f j (t).

B. LOCAL CONTROLLERS FOR SUB-SYSTEMS WITH FLOW CONVERGENT TOPOLOGY 1) Dynamic model of sub-systems with flow convergent topology
As stated in Definition 2, sub-systems with flow convergent topology have n source reservoirs and one receptor reservoir (see Figure 2(b)). For convention, we use the indices 1, . . . , n, for identifying source reservoirs, while the receptor reservoir is marked with the index n + 1.
The differential equations describing the change along the time of the source reservoirs' wastewater volume arė where k i is a parameter that depends on the geometry of the i-th source reservoir, and x i (t) ∈ [0, 1] is the opening percentage of the valve that connects the i-th source reservoir with the receptor reservoir.
On the other hand, the dynamic model of the receptor reservoir's wastewater volume is given bẏ external outflow (9) Figure 7 shows a schematic representation of a sub-system with flow convergent topology. Notice that all the involved variables and parameters of the dynamic model described before are depicted in this figure.
v n+1 (t) q in n (t) Figure 7. A sub-system with convergent topology.
Remark 2. The variables of the divergent topology model (see Figure 6) are the same as the ones used in the divergent case (see Figure 7). Nonetheless, the equations that describe both models ( (4) and (5) for divergent topologies, and (8) and VOLUME 4, 2016 (9) for convergent ones) are different. This fact is because while in the convergent model we have n receptor reservoirs and one source, in the divergent model we have n sources and only one receptor reservoir.

2) Replicator dynamics-based controller for flow convergent topology
Let us define the parameters of the replicator dynamics-based controller for each sub-system having a flow convergent topology.
Similarly as in the flow divergent case, we have a strategy set given by S = {1, . . . , n + 1}. The population playing the i-th strategy, where i = 1, . . . , n, is the opening percentage of the i-th source tank's valve, i.e., x i (t). Besides, the population playing the strategy n + 1, i.e., x n+1 (t), is equal to the sum of the remaining opening percentages of all valves. Therefore, the population mass is again equal to n.
Fitness functions must capture the fact that if the remaining volume of a source tank is low, then it is more appealing to increase the opening percentage of its valve. Having this idea in mind, we define the fitness functions as minus the remaining volume of each tank, i.e., As in the flow divergent case, we have added the same barrier term εb i (x i (t), 1) to the first n fitness functions.
Finally, the replicator dynamics-based controller for subsystems with flow convergent topology is given by the equations in (7), with the same average fitness.
Remark 3. The local controller of each sub-system handles the valves that belong to the subsystem. Therefore, if a valve (or what is the same, an edge of the graph that models the UDS) belongs to more than one sub-system, there would be ambiguity because more than one controller would make decisions about the valve opening percentage. This is the reason for the partitioning constraint C3 (see Section IV).

VI. STABILITY ANALYSIS
This section analyzes the stability of each sub-system controlled via replicator dynamics. We focus our analysis in a special case where the opening percentages of the subsystems' valves satisfy the following assumptions: . , x * n , be the opening percentages of a sub-system's valves at equilibrium. The value x * i > 0, for all i = 1, . . . , n. Thus, none of the sub-system's valves are fully closed at equilibrium.

Assumption 3. The value of x *
i is such that the barrier terms in (6) and (10) can be neglected in the computation of the fitness functions. Hence, x * i is far from being fully opened. Furthermore, we make the following assumption on the external inflows and outflows of sub-systems. Assumption 4. Consider a sub-system having either convergent or divergent topology. External inflows and coefficients related to external outflows are constants, i.e., q in i (t) =q in i and γ out i (t) =γ out i , for all i = 1, . . . , n + 1.

1) Equilibrium point
Under Assumptions 2-4, the equilibrium point of a subsystem with flow divergent topology controlled via population dynamics has an important property, which is given in the following proposition.

Proposition 1. Let Assumptions 2-4 hold. If
n+1 i=1 x i (0) = n, then the equilibrium point of the closed-loop system given by (4), (5) and (7) is where v * i is the wastewater volume of the i-th reservoir at equilibrium. Notice that the remaining capacities of all reservoirs are the same at equilibrium, i.e.,v , for all i = 1, . . . , n + 1.
Proof. Equilibrium condition for (4) and (5) is given by On the other hand, due to Assumption 2, all fitness functions reach the same value at equilibrium (see (7)). Moreover, fitness functions are equal to the remaining capacity of each reservoir for Assumption 3 (see (6)). Therefore, all subsystem's reservoirs have the same remaining capacity at equilibrium, i.e.,v i −v * i =v j −v * j , for all i, j = 1, . . . , n+1. Using this fact and after some algebraic manipulations on (12), we obtain the first equations of the equilibrium point given in (11), i.e., v * i , for all i = 1, . . . , n + 1, and x * i , for all i = 1, . . . , n. The last equation of (11), i.e., the equilibrium state x * n+1 , is computed using the fact that The fact that all reservoirs have the same remaining capacity at equilibrium implies a fairly distribution of wastewater throughout the sub-system. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3177631, IEEE Access

2) Stability
The next result guarantees the stability of the equilibrium point given in Proposition 1 for a sub-system with divergent topology controlled via replicator dynamics.
Proof. The sub-system controlled via replicator dynamics, which is given in (4), (5) and (7), can be written in error coordinates with respect to the equilibrium point in (11) as follows: Σ2 : where e vi (t) = v i (t) − v * i and e xi (t) = x i (t) − x * i , for all i = 1, . . . , n + 1. Notice that Σ 1 is the UDS sub-system and Σ 2 is the controller.
We use Lyapunov theory to prove stability of the equilibrium point. For Σ 1 , consider the positive definite function where e v (t) = e v1 (t), . . . , e vn+1 (t) ⊤ . The derivative of For computingV 1 e v (t) , we have used the fact that n+1 i=1 x i (t) = n, for all t ≥ 0 (Lemma 1 guarantees this property since n+1 i=1 x i (0) = n by assumption). The function Ψ 1 e v (t) is negative definite since x i (t) ∈ [0, 1] (this constraint also follows from Lemma 1 and by the effect of the barrier functions given in (6)) and k n+1 < 4γ out i by assumption.
On the other hand, consider the following positive definite function (which is based on the relative entropy function [53]) for Σ 2 : where e x (t) = e x1 (t), . . . , e xn+1 (t) ⊤ . Positive definiteness of V 2 e x (t) can be proven using Jensen's inequality as follows: since the logarithm function is strictly concave and given the fact that In addition, according to Jensen's inequality, V 2 e x (t) = 0 only at x i (t) = x * i , for all i = 1, . . . , n + 1. Therefore, V 2 is positive definite.
The derivative of V 2 e x (t) along the trajectories of Σ 2 iṡ The function Ψ 2 e x (t) ≤ 0 since b i (x i (t), 1) is monotonically decreasing with respect to x i (t). Furthermore, although the value of Ψ 2 e x (t) does not depend on e xn+1 (t), only if e xi (t) = 0, for all i = 1, . . . , n. This condition implies that e xn+1 (t) = 0 (which follows from the fact that e xn+1 (t) = − n i=1 e xi (t) according to Lemma 1). Based on the results described above, we propose the Lyapunov function candidate V e v (t), e x (t) = V 1 e v (t) + Ψ 2 e x (t) , which is negative definite. Thus, we can conclude that the equilibrium point in (11) is asymptotically stable.
Theorem 1 is relevant since it guarantees an appropriate performance of the proposed controller. Notice that the result stated in this theorem implies that if a reservoir has a lower remaining capacity than the other reservoirs (e.g., if a reservoir is going to overflow), then the controller tries to mitigate this undesirable situation by reallocating the wastewater in all reservoirs of the sub-system. This wastewater reallocation is performed in such a way that all reservoirs have the same remaining capacity.

B. FLOW CONVERGENT TOPOLOGY CASE
For sub-systems with flow convergent topology, we have similar results than those shown for the flow divergent case.

1) Equilibrium point
As stated in the next proposition, a sub-system with flow convergent topology controlled via replicator dynamics has an equilibrium point in which the remaining capacities of all reservoirs are the same.
(15) Notice that, in this equilibrium point, the remaining capacities of all reservoirs (i.e.,v i − v * i , for all i = 1, . . . , n + 1) are the same.
Proof. The proof follows similar ideas that the ones used in the proof of Proposition 1.
As in the flow divergent case, the equilibrium point in (15) implies a well-balanced wastewater distribution among the reservoirs of the sub-system. Indeed, notice that if all reservoirs have the same remaining capacity, there are no overflows in any reservoir.

2) Stability
The next theorem shows that the replicator dynamics based controller drives the sub-system to the equilibrium point given in (15). Theorem 2. Let Assumptions 2-4 hold. In addition, letv and k be positive real values, and assume thatv i =v, for all i = 1, . . . , n + 1, and k i = k, for all i = 1, . . . , n. Furthermore, assume that k < 4γ out n+1 n . If the initial conditions of the controller in (7) satisfy n+1 i=1 x i (0) = n and x i (0) > 0, for all i = 1, . . . , n + 1, then the equilibrium point stated in (15) is asymptotically stable under the dynamics of the closedloop system given in (7)- (5).
Proof. Since all reservoir volumes are the same by assumption, then there exists a positive real value v * such that v * i = v * , for all i = 1, . . . , n + 1 (see (15)). Taking where V 2 e x (t) is defined in (14), as Lyapunov function candidate and following the same procedure of the proof of Theorem 1, we can conclude that the derivative of V e v (t) along the trajectories of the closed-loop system is negative definite. Therefore the equilibrium point is asymptotically stable.
For the case of flow convergent topology, notice that we impose additional constraints (compared to the flow divergent case) to guarantee stability of the desired equilibrium point. Specifically, we require that all reservoirs have the same volumesv i and the same coefficients k i . These constraints are met for sub-systems where all reservoirs have the same size and geometry.
Once the stability has been proven for both divergent and convergent topology sub-systems, we test the proposed control strategy in a case study.

VII. CASE STUDIES AND DISCUSSION
We test the performance of the proposed decentralized controller in two simulation scenarios. Furthermore, we compare our method with two techniques that have been employed in real-time control for UDSs: decentralized LQR with constraints (e.g., see [14], [54]) and centralized MPC (e.g., see [55]).
For synthesizing the decentralized LQR scheme, we use the partitioning algorithm described in Section IV and design local LQRs per each sub-system. On the other hand, for designing the MPC strategy, we use the information of the whole UDS to compute the control signals. Hence, we notice that the technique developed in this article and the decentralized LQR scheme have disadvantages compared to MPC in terms of the information availability (i.e., local vs global). However, there are two reasons for using centralized MPC in our comparison. First, we want to compare our method with another that might guarantee an optimal performance as is the case of MPC. Second, the global information requirement of centralized MPC can be relaxed by decoupling the optimization problem that MPC solves at each step and then using a distributed optimization algorithm (e.g., ADMM [43]). This relaxation allows MPC to use local information at the cost of increasing the computation time since distributed optimization algorithms have slower convergence rates than their centralized counterparts.

A. CASE STUDY 1
The first scenario corresponds to the UDS shown in Figure  4(a). This UDS has seven reservoirs and four sub-systems (as explained in Section IV-D). Half of the sub-systems has flow divergent topology and the other half has flow convergent topology (see Figure 5). For this case study, all reservoirs have the same parameters, i.e., k i = 0.005 s −1 and v i = 274.80 m 3 , for all i = 1, . . . , 7 (these parameters have been taken from [23]). We simulate a rainfall scenario where flooding typically occurs. The direct runoff hydrograph for each reservoir, which models the external inflows, is shown in Figure 8(c). These runoff hydrographs are Gaussian functions whose parameters are randomly chosen as follows: center is in the interval 0.5-0.75 hours, width is in the interval 0-0.5 hours, and amplitude is in the interval 0-0.6 m 3 /s. We also assume that there are no leakages at reservoirs. Therefore, external outflows of reservoirs 1, 2, 4, 5, and 6 are only due to wastewater flowing to other sub-systems of the UDS, while external outflows of reservoirs 3 and 7 correspond to the wastewater flowing outside the UDS (e.g., to wastewater treatment plants or to water bodies). In this case, we assume that γ out 3 (t) = γ out 7 (t) = 0.005 s −1 . Evolution of the remaining capacities of each UDS reservoir is shown in Figure 8(a). In this case, all valves of the system are completely opened without any controller acting on them. If the remaining capacity of some reservoir reaches the dotted line, it indicates that the reservoir is full of wastewater. On the other hand, if the remaining capacity of some reservoir is below the dotted line, it implies that wastewater exceeds the capacity of the reservoir, and thus flooding occurs. Notice that, for the uncontrolled case shown in Figure 8(a), there is an overflow in reservoir 7. In fact, the total overflow is around 71.3 m 3 . This undesirable situation is avoided by using the decentralized controllers based on replicator dynamics with parameters β = 0.0006, ε = 0.0001, and barrier functions b i = 1 xi(t)−1 . Evolution of the remaining capacities of reservoirs under the proposed control strategy is shown in Figure 8(b), while the control signals (which correspond to the opening percentages of system's valves) are depicted in Figure 8(d). Notice that, under the replicator dynamics-based controller, there are no overflows in the UDS. The reason for this improved performance is that the wastewater is properly allocated among all the UDS reservoirs. The key point of the achieved wastewater allocation is that the proposed controller takes advantage of the capacities of upstream reservoirs. In the uncontrolled scenario, upstream reservoirs are promptly emptied causing downstream reservoirs receive a large amount of wastewater. Hence, reservoir 7 gets overloaded. On the contrary, when the controller is applied, valves of upstream reservoirs are closed (e.g., see valve 1 and valve 3 in Figure 8(d)). Therefore, these reservoirs are more filled than in the uncontrolled scenario causing that the flow downstream is relieved and flooding is avoided.
For this case study, we also test a decentralized LQR scheme and a centralized MPC synthesis. For both controllers the cost function is defined by i v ⊤ i v i , where the sum is over the reservoirs that belong to the corresponding subsystem in the case of LQR, and over all the reservoirs of the UDS in the case of MPC. Notice that this cost function penalizes fuller reservoirs. For both controllers, we constraint the opening percentage of valves to belong to the interval [0, 1], i.e., x i ∈ [0, 1]. Moreover, the sample time is 5 seconds for LQR, while for MPC, the sample time is 1 minute, the prediction horizon is 60 steps, and the control horizon is 10 steps. The performance of the controlled UDS is shown in Figure 9. Capacities of the UDS's reservoirs are depicted in Figure 9(a) (for LQR) and in Figure 9(b) (for MPC). In addition, control signals are shown in Figure 9(c) (for LQR) and in Figure 9(d) (for MPC).
Using the decentralized LQR, flooding occurs. However, it is significantly lower than the one obtained without using any controller. Specifically, the flooding is 18.2 m 3 , which corresponds to a decrease of flooding of 75% compared to the baseline scenario (no control). Although the system behaves better under LQR than without control, the performance of LQR is worse than the decentralized controller based on replicator dynamics, which completely eliminates flooding under the same simulation conditions. MPC also avoids flooding. Indeed, it is worth noting that the control signals generated when the centralized MPC is employed (Figure 9(d)) are quite similar to those produced under the decentralized controller based on replicator dynamics. This fact occurs even when the centralized MPC uses global information (i.e., information of the whole UDS) to compute the control inputs while the replicator dynamics only employ local information (i.e., information of each subsystem).
Furthermore, we want to highlight that both LQR and MPC strategies try to equalize the remaining capacities of reservoirs (see Figures 9(a) and 9(b)) in the same way as the replicator dynamics do. This happens even though neither LQR nor MPC controllers are designed using a cost function that explicitly penalizes the differences between the remaining capacities of reservoirs (the employed cost function only penalizes fuller reservoirs as was explained when the design parameters of LQR and MPC were presented).
Finally, to provide a quantitative comparison of the reduction of flooding using the different control techniques, i.e., local controllers based on replicator dynamics, decentralized LQR and centralized MPC, we simulate the UDS of case study 1 under the same conditions described before, but increasing the rain intensity. To do so, we amplify the peak of the Gaussian functions that model the rain by the factors given in Table 1. The obtained flooding is also shown in that table. Notice that comparing both decentralized strategies, the controller based on replicator dynamics outperforms LQR. For instance, in the heaviest rain scenario, the percentage of flooding reduction under replicator dynamics is 51% compared to LQR. On the other hand, although centralized MPC shows the best performance, the decentralized control strategy based on replicator dynamics keeps flooding at levels VOLUME  (c) (d) close to the optimal. At this point, it is worth mentioning that the replicator dynamics based controller benefits from lower computation requirements compared to MPC. In fact, the replicator dynamics use a simple model given by a set of firstorder ordinary differential equations, while MPC requires solving an optimization problem at each time step. Moreover, features such as modularity provided by a decentralized scheme also represent a key advantage from the controller based on replicator dynamics with respect to the centralized MPC approach. Finally, the need of merging information of the entire system in a centralized controller could represent a drawback with respect to a non-centralized control approach, which does not need such a strong condition over the information availability.

B. CASE STUDY 2
We also analyze the performance of the proposed control technique in a case study based on a real scenario. To do that, we simulate a portion of the Bogotá (Colombia) stormwater UDS whose description is presented in [37]. Figure 10 shows the structure of this UDS, which has 16 reservoirs. Parameters of each reservoir are given in Table 2.
Simulation conditions reproduce a heavy rains scenario, where the external inflow of each reservoir is due to direct runoff. Figure 11 shows the hydrographs that models these runoffs. Similarly as in the previous case study, hydrographs are Gaussian functions whose parameters lie in the following intervals: center is between 0.5-0.75 hours; width is between  (c) (d) 0-0.5 hours; and amplitude is between 0-0.22 m 3 /s. Figure  12 shows the behavior of the UDS without control (i.e., reservoirs' valves are completely open during all the time) for a time window of 3.5 hours. Throughout the simulation, wastewater only fills half of most reservoirs. In contrast, reservoirs 12 and 14 are overloaded during the first half of the simulation. Indeed, the remaining capacity of reservoir 14 is negative during this first half implying that flooding occurs. The total flooding is around 797.1 m 3 . This situation is mainly due to two facts: first, reservoir 14 is a downstream reservoir, which receive wastewater from a big portion of the UDS; second, reservoir 14 is the one with the lowest coefficient k i , i.e., it has the lowest capacity to evacuate wastewater (according to (8), the larger the coefficient k i , the larger the outflow of the i-th reservoir trough its output valve).
Once the baseline (corresponding to the behavior of the system without control) has been established, we implement the control strategy proposed in Section V. First, we apply the partitioning algorithm of Section IV and identify 6 sub-systems with flow convergent topology (each sub-system is associated with a color in Figure 10). Table 3 summarizes the structure of the identified sub-systems. For each one of them, we implement a replicator dynamics based controller with parameters β = 0.0003, ε = 5, and barrier functions 1 xi(t)−1 . Simulation results of the controlled UDS are shown in Figure 13. Additionally, the control inputs (opening percentage of the reservoirs' valves) are shown in Figure 14.
The first thing that we notice in Figure 13 is that, unlike the uncontrolled case, there are no big differences in the remaining capacities of reservoirs. This behavior is due to the fact that the replicator dynamics based controllers seek to equalize the remaining capacities as proven in Section VI. For instance, see the behavior of reservoirs 14 (brown line) and 15 (orange line), and their corresponding valves. From hour 1-2, reservoir 14 has lower remaining capacity than reservoir 15. Thus, during this time, the valve that allows reservoir 14 to evacuate wastewater is almost fully opened while valve of reservoir 15 is closed. This fact allows reservoirs 14 and 15 to reach almost the same remaining  Rain at each reservoir capacity from hour 2 onwards. Summarizing, when using the local controllers based on replicator dynamics, the total capacity of the UDS is better used because wastewater is stored in underloaded reservoirs. On the other hand, the total   14, 15 16 As in the first case study, we implement a decentralized LQR and a centralized MPC to compare the performance of the controller based on replicator dynamics. The design process of LQR and MPC is the same described in case study 1. Besides, the sample time of LQR is 2 seconds, while MPC parameters are chosen as follows: sample time is 10 seconds, prediction horizon is 360, and control horizon is 60. Simulation results under LQR are shown in Figure 15 (reservoirs' capacities) and Figure 16 (control signals). In addition, simulation results under MPC are shown in Figure  17 (reservoirs' capacities) and Figure 18 ( Control inputs using RD-based controllers performed in the first case study, the results are summarized in Table 4. Notice that we obtain similar results as the ones of case study 1, i.e., replicator dynamics controller outperforms LQR and its behavior is quite close to centralized MPC, but with computational advantages as discussed next.

C. COMPUTATIONAL BURDEN COMPARISON
An important advantage of the proposed control technique is its low computational burden, especially when we compare replicator dynamics with optimization-based controllers such as LQR and MPC. This advantage rests on the philosophy behind the controllers. While techniques based on optimization should minimize a cost function subject to constraints at each control step, replicator dynamics only require the solution of a set of ordinary differential equations. The key point is that iterative algorithms are required for solving constrained optimization problems. This fact implies that, for each control step, optimization-based controllers should perform several iterations, which increases the computational burden. In contrast, each time step of replicator dynamics only requires the evaluation of the right-hand side of the replicator equation in (7) to update the control input (using, for instance, the Euler approximation). Clearly, the computational complexity of replicator dynamics is significantly lower. To demonstrate this fact, we make a comparison of the computational time taken for simulating each case study under the three different control strategies implemented in Sections VII-A and VII-B, i.e., decentralized replicator dynamics, decentralized LQR and centralized MPC. To perform the simulations, we use a machine whose characteristics are summarized in Table 5.
The results of the comparison are shown in Table 6. In this table, the sub-index of r indicates the rain scale factor, which is explained in the simulation conditions of both case studies (see Sections VII-A and VII-B). For instance, r 1.10 denotes that the rain amplitude of the baseline case is increased 10%.
Notice that in all cases the controller based on replicator dynamics is the fastest. In fact, for case study 1, replicator dynamics is three times faster than decentralized LQR and 14 times faster than centralized MPC (the slower method). These differences are higher for case study 2, where replicator dynamics is six times faster than decentralized LQR and 122 times faster than centralized MPC. On the other  [56] hand, it is worth mentioning that replicator dynamics exhibit better scalability. Indeed, if we study the computational time taken in case study 1 (UDS with 7 reservoirs) compared to case study 2 (UDS with 16 reservoirs), it can be noticed that the increment is 32 times for centralized MPC, seven times for decentralized LQR and only 3.7 times for replicator dynamics. The fact that LQR and replicator dynamics scale much better than MPC is a consequence of the decentralized nature of the first two controllers, whose parallelization plays a key role (i.e., each UDS sub-system is capable of computing its own control inputs independently). In addition, another factor that makes replicator dynamics beat the two optimization-based methods in terms of scalability is the different paradigm that replicator dynamics use to update the control inputs. Specifically, when the UDS scale increases, optimization-based controllers require to solve minimization problems with a larger number of decision variables (under this situation, centralized MPC is highly penalized because it must deal with all decision variables of the whole UDS). On the contrary, replicator dynamics only need the evaluation of more algebraic equations to update the control input. As explained before, this evaluation requires significantly lower computations than solving a constrained optimization problem.

VIII. CONCLUSIONS
Floods can be mitigated using real time control of Urban Drainage Systems (UDSs). Since the size of these kinds of systems is generally of large-scale nature, designing a central control (which drives all the actuators and has information from all UDS's sensors) is problematic. Instead, it is desirable to design decentralized control schemes, i.e., schemes with multiple controllers, where each one of them has only partial information of the system. We have shown that graph models and population dynamics can be employed to design local controllers that independently handle partitions of the UDS. The proposed scheme allocates wastewater among UDS collectors (modeled as reservoirs) in an efficient way. Specifically, we have proven that the designed local controllers perform wastewater allocation in such a way that the collectors of the UDS tend to reach the same remaining storing capacity while preserving closed-loop system stability. Remaining capacities of reservoirs using decentralized LQR Figure 15. Evolution of the remaining capacities of reservoirs of the Bogotá (Colombia) case study under decentralized LQR.    lations have shown that the performed wastewater allocation mitigates flooding even in the face of heavy rain scenarios.
As future work, we propose to improve the performance of the controlled system by incorporating prediction of the UDS state and forecasting of external inputs. Besides, we propose to study tuning algorithms for finding the optimal parameters of the local controllers.