Allocation is the process of assigning an allocated flow rate to demand nodes in the physical layer of the model based on information about sources, the different demand nodes over various priorities, constraints introduced by nodes, local water availability and graph topology. The allocation procedure implemented in Ribasim is heavily inspired by the maximum flow problem.
The allocation problem is solved per subnetwork (and main network) of the Ribasim model. Each subnetwork is used to formulate an optimization problem with the JuMP package, which is solved using the HiGHS solver. For more in-depth information see also the example of solving the maximum flow problem with JuMP.jlhere.
Before the optimization for each priority there is a simple step that tries to allocate flow to the UserDemand nodes from the Basin directly upstream.
Note
within this Allocation section the main network is also considered to be a subnetwork.
2 The high level algorithm
The allocation algorithm contains 3 types of optimization:
internal_sources, where flows are allocated within a subnetwork by only using sources inside the subnetwork;
collect_demands, where flows are allocated within a subnetwork by only using the main network inlet(s) as a source, with demands reduced by allocations in internal_sources. The allocated flows in this optimization type are not used. The goal is to see the flow through the main network inlet(s), which is interpreted as the subnetwork demand;
allocate, where all available sources are used and the final allocated flows for the users are determined.
The full algorithm goes through the following steps:
Perform internal_sources followed by collect_demands for all subnetworks apart from the main network;
Perform allocate for the main network;
Perform allocate for the other subnetworks.
If no main network is present, then 1 and 2 are skipped.
3 Elements of allocation
The following data of the parameters and state of a Ribasim model are relevant for the allocation problem.
3.1 Schematisation input
3.1.1 The subnetwork
The allocation problem is solved per subnetwork, which is given by a subset \(S \subset V\) of node IDs. Different subnetworks are disjoint from eachother. Nodes can also not be part of any subnetwork.
3.1.2 Source flows
Sources are indicated by a set of edges in the subnetwork \[
E_S^\text{source} \subset E,
\] which are automatically inferred as all edges that point out of LevelBoundary or FlowBoundary nodes. That is, if \((i,j) \in E_S^\text{source}\), then the average over the last allocation interval \(\Delta t_{\text{alloc}}\) of the of the flow over this edge \[
\frac{1}{\Delta t_{\text{alloc}}}\int_{t - \Delta t_{\text{alloc}}}^tQ_{ij}(t') dt'
\] is treated as a source flow in the allocation problem. These edges are either coming from a boundary/source node (e.g. a level or flow boundary) or connect the main network to a subnetwork. For the definition of \(Q_{ij}\) see the formal model description.
3.1.3 User demands
The subnetwork contains a subset of UserDemand nodes \(U_S \subset S\), who all have static or time varying demands over various priorities \(p\): \[
d^p_i(t), \quad i \in U_S, p = 1,2,\ldots, p_{\max}.
\]
Note
On this page we assume that the priorities are given by all integers from \(1\) to some \(p_{\max} \in \mathbb{N}\). For the Ribasim input this is not a requirement; some of these in between priority values can be missing, only the ordering of the given priorities is taken into account.
3.1.4 Flow demands
The subnetwork contains a subset of nodes \(FD_S \subset S\) which have a demand of a single priority \(p_{\text{fd}}\) for the flow through that node. With this we define \[
d^p_i(t) =
\begin{cases}
0 \text{ if } p \ne p_{\text{fd}} \\
d^{p_{\text{df}}} \text{ if } p = p_{\text{fd}}
\end{cases}
\] for all \(i \in FD_S\). Here \(d^{p_{\text{df}}}\) is given by the original flow demand minus the flows trough node \(i\) at all priorities \(p < p_{\text{fd}}\).
3.2 Simulation (physical layer) input
3.2.1 Vertical fluxes and local storage
Apart from the source flows denoted by edges, there are other sources of water in the subnetwork, associated with the Basins in the subnetwork \(B_S = B \cap S\). First, there is the average over the last allocation interval \(\Delta t_{\text{alloc}}\) of the vertical fluxes (precipitation, evaporation, infiltration and drainage) for each Basin: \[
\phi_i(t) = \frac{1}{\Delta t_{\text{alloc}}}\int_{t - \Delta t_{\text{alloc}}}^t \left[Q_{P,i}(t') - Q_{E,i}(t') + Q_{\text{drn},i}(t') - Q_{\text{inf},i}(t') \right] dt', \quad \forall i \in B_S.
\]
We consider fluxes into the basin to be positive and out of the basin to be negative. For more information see the natural water balance terms.
Secondly, there is either a supply or demand from the storage in the basin. Given a minimum level \(\ell_{\min, i}\) and a maximum level \(\ell_{\max, i}\) which correspond to a minimum storage \(s_{\min, i}\) and maximum storage \(s_{\max, i}\) respectively, we get a flow supply of \[
F^{\text{basin out}}_{\max, i} = \max\left(0.0, \frac{u_i(t)-s_{\max,i}}{\Delta t_{\text{alloc}}} + \phi_i(t)\right)
\]
and a demand of \[
d^p_i = \max\left(0.0, \frac{s_{\min,i} - u_i(t)}{\Delta t_{\text{alloc}}} - \phi_i(t)\right),
\]
for all \(i \in B_S\). Note that the basin demand has only a single priority, so for other priorities this demand is \(0\).
3.2.2 Constraining factors
3.2.2.1 Flow magnitude and direction constraints
Nodes in the Ribasim model that have a max_flow_rate, i.e. Pump, Outlet and LinearResistance, put a constraint on the flow through that node. Some nodes only allow flow in one direction, like Pump, Outlet and TabulatedRatingCurve.
3.2.2.2 UserDemand return flows
UserDemand nodes dictate proportional relationships between flows over edges in the subnetwork. The return factor is given by \(0 \le r_i(t) \le 1, i \in U_S\).
3.3 The subnetwork
The subnetwork consists of a set of nodes \(S \subset V\) and edges
\[
E_S = (S \times S) \cup E_S^\text{source},
\]
i.e. the edges that lie within the subnetwork together with the source edges (which can be partially outside the subnetwork). The nodes in \(S\) together with the connected nodes outside the subnetwork are called the extended subnetwork.
3.3.1 Capacities
Each edge in the subnetwork has an associated capacity. These capacities are collected in the sparse capacity matrix \(C_S \in \overline{\mathbb{R}}_{\ge 0}^{n\times n}\) where \(n\) is the number of nodes in the extended subnetwork. An edge capacity is infinite if there is nothing in the model constraining the capacity.
The capacities are determined in different ways:
If an edge does not exist in the allocation network, i.e. \((i,j) \notin E_S\) for certain \(1 \le i,j\le n'\), then \((C_S)_{i,j} = 0\);
The capacity of the edge \(e \in E_S\) is given by the smallest max_flow_rate of the nodes along the equivalent edges in the subnetwork. If there are no nodes with a max_flow_rate, the edge capacity is infinite;
If the edge is a source, the capacity of the edge is given by the flow rate of that source;
There are also capacities for special edges:
\(C^{LD}_S \in \mathbb{R}^b_{\ge 0}\) where \(b = \# B_S\) is the number of basins, for the flow supplied by basins based on level demand (this capacity is 0 for basins that have no level demand).
\(C^{FD}_S \in \mathbb{R}^c_{\ge 0}\) where \(c = \# FD_S\) is the number of nodes with a flow demand, for the flow supplied by flow buffers at these nodes with a flow demand.
\(C^{UD}_S \in \mathbb{R}^f_{\ge 0}\) where \(f = \# U_S\), for the flow supplied by the user demand outflow source whose capacity is given by return flows.
4 The optimization problem
The optimization problem for a subnetwork consists of a quadratic objective function with associated linear constraints on a set of variables, all of which are introduced below.
4.1 The optimization variables
There are several types of variable whose value has to be determined to solve the allocation problem:
The flows \(F \in \mathbb{R}_{\ge 0}^{n\times n}\) over the edges in the allocation network;
The flows \(F^\text{basin out}_{i}, F^\text{basin in}_{i} \geq 0\) for all \(i \in B_S\) supplied and consumed by the basins with a level demand respectively;
The flows \(F^\text{buffer out}_{i}, F^\text{buffer in}_{i} \ge 0\) for all \(i \in FD_S \cup FF_S\) supplied and consumed by the flow buffers of nodes with a flow demand.
4.2 The optimization objective
The goal of allocation is to get the flow to nodes with demands as close as possible to these demands. To achieve this, a sum of error terms is minimized.
The error between the flows and user demands is denoted by \(E_{\text{user demand}}\), where \[
E_{\text{user demand}} = \sum_{(i,j)\in E_S\;:\; i\in U_S} d_j^p(t)\left(1 - \frac{F_{ij}}{d_j^p(t)}\right)^2
\]
Note
When performing main network allocation, the connections to subnetworks are also interpreted as UserDemand nodes with demands determined by subnetwork demand collection.
This type of objective cares about the fraction of the demand allocated, and will lead to an equal fraction of all demands allocated when possible. For a discussion on this see here.
Likewise, the error of level demands from basins is the squared relative difference between flows consumed by basins and basin demands. \[
E_{\text{level demand}} = \sum_{i \in B_S} d_i^p(t)\left(1 - \frac{F_i^\text{basin in}}{d_i^p(t)}\right)^2
\]
Lastly, the error of the flow demands is given as below. \[
E_{\text{flow demand}} = \sum_{i \in FD_S} d_i^p(t)\left(1 - \frac{F_i^\text{buffer in}}{d_i^p(t)}\right)^2
\]
for the set of in-neighbors and out-neighbors of a node in the network respectively.
Flow conservation: For all nodes \(k\) that are not a source or a sink (i.e. FlowBoundary, LevelBoundary, UserDemand) we have a flow conservation constraint: \[
\sum F_{\text{out special}} + \sum_{j \in V^{\text{out}}_S(k)} F_{kj} = \sum F_{\text{in special}} + \sum_{i \in V^{\text{in}}_S(k)} F_{ik}, \quad \forall k \in B_S.
\tag{1}\]
In here, we have the following special flows:
If \(k\) is a basin with a flow demand, there is a special outflow \(F^{\text{basin in}}_k\) and a special inflow \(F^{\text{basin out}}_k\);
If the node has a buffer (see here) there is a special outflow \(F^{\text{buffer in}}_k\) and a special inflow \(F^{\text{buffer out}}_k\).
Note
In the above, the placement of the basin and buffer flows might seem counter-intuitive. Think of the storage or buffer as a separate node connected to the node with the demand.
Capacity: the flows over the edges are bounded by the edge capacity: \[
F_{ij} \le \left(C_S\right)_{ij}, \quad \forall(i,j) \in E_S.
\tag{2}\] By the definition of \(C_S\) this also includes the source flows. The same holds for the basin outflows:
There are several things to note about the source constraints: - The sources are not all used at once. There is an optimization for each source in a subnetwork, where only one source has nonzero capacity. - When performing subnetwork demand collection, these capacities are set to \(\infty\) for edges which connect the main network to a subnetwork.
Similar constraints hold for the flow out of basins, flow demand buffers and user demand outflow sources: \[
F^\text{basin out}_{i} \le (C^{FD}_S)_i, \quad \forall i \in B_S,
\]
\[
F_{ij} \le (C^{UD}_S)_i, \quad \forall i \in U_S, \quad V_S^{\text{out}}(i) = \{j\}.
\] Here we use that each UserDemand node in the allocation network has a unique outflow edge. The user outflow source capacities are increased after each optimization solve by the return fraction: \[
r_i(t) \cdot F_{ki}, \quad V_S^{\text{in}}(i) = \{k\}.
\]
Flow sign: Furthermore there are the non-negativity constraints for the flows and allocations, see The optimization variables.
4.4 Example
The following is an example of an optimization problem for the example shown here:
Code
usingRibasimusingRibasim: NodeIDusingSQLiteusingComponentArrays: ComponentVectortoml_path =normpath(@__DIR__, "../../generated_testmodels/allocation_example/ribasim.toml")model = Ribasim.Model(toml_path)(; u, p, t) = model.integratorallocation_model = p.allocation.allocation_models[1]priority_idx =1Ribasim.set_initial_values!(allocation_model, u, p, t)Ribasim.set_objective_priority!(allocation_model, u, p, t, priority_idx)println(p.allocation.allocation_models[1].problem)