Allocation

1 Introduction

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.jl here.

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:

1. Perform internal_sources followed by collect_demands for all subnetworks apart from the main network;
2. Perform allocate for the main network;
3. 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.

3.1.2 Source flows

Sources are indicated by a set of edges in the subnetwork \[ E_S^\text{source} \subset E. \] 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}}\). 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\). Firstly 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 FractionalFlow and UserDemand return flows

Both FractionalFlow and UserDemand nodes dictate proportional relationships between flows over edges in the subnetwork. UserDemands have a return factor \(0 \le r_i \le 1, i \in U_S\), and FractionalFlow nodes have an associated fraction \(f_j\) which can be updated by DiscreteControl.

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;
• If an edge comes from a node with a flow demand, it has infinite capacity at priorities other than this of this flow demand, and zero capacity otherwise.

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 is a linear optimization problem consisting of an objective function with associated 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 or fractional flow outneighbors.

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 error of terms is minimized.

\[ \min E_{\text{user demand}} + E_{\text{level demand}} + E_{\text{flow demand}} \]

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} \left| F_{ij} - d_j^p(t)\right| \]

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 absolute amount of water allocated to a demand. It treats all deviations equally which means it doesn’t give larger punishment per flow unit if deviations increase.

The absolute value applied here is not supported in a linear programming context directly; this requires introduction of new variables and constraints. For more details see here.

Likewise, the error of level demands from basins is the absolute difference between flows consumed by basins and basin demands. \[ E_{\text{level demand}} = \sum_{i \in B_S} \left| F_i^\text{basin in} - d_i^p(t)\right| \]

Lastly, the error of the flow demands is given as below. \[ E_{\text{flow demand}} = \sum_{i \in FD_S} \left| F_i^\text{buffer in} - d_i^p(t)\right| \]

4.3 The optimization constraints

For convenience, we use the notation

\[\begin{align} V^{\text{out}}_S(i) = \left\{j \in V : (i,j) \in E_S\right\} \\ V^{\text{in}}_S(j) = \left\{i \in V : (i,j) \in E_S\right\} \end{align}\]

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:

\[ F^{\text{basin out}}_{i} \le F^{\text{basin out}}_{\max, i}, \quad \forall i \in B_S. \]

Note

When performing subnetwork demand collection, these capacities are set to \(\infty\) for edges which connect the main network to a subnetwork. For all other sources the capacity is set to \(0\), so that demand collection only uses flow from the main network inlet.

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^\text{buffer out}_{i} \le (C^{FD}_S)_i, \quad \forall i \in FD_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 \cdot F_{ki}, \quad V_S^{\text{in}}(i) = \{k\}. \]

• Fractional flow: Let \(L_S \subset V_S\) be the set of nodes in the max flow graph with fractional flow outneighbors, and \(f_j\) the flow fraction associated with fractional flow node \(j \in V_S\). Then \[ F_{ij} \le f_j \sum_{k\in V^\text{in}_S(i)} F_{ki} \qquad \forall i \in L_S, \; j \in V_S^\text{out}(i). \tag{3}\]

• 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:

using Ribasim
using Ribasim: NodeID
using SQLite
using ComponentArrays: ComponentVector

toml_path = normpath(@__DIR__, "../../generated_testmodels/allocation_example/ribasim.toml")
p = Ribasim.Model(toml_path).integrator.p
u = ComponentVector(; storage = zeros(length(p.basin.node_id)))

allocation_model = p.allocation.allocation_models[1]
t = 0.0
priority_idx = 1

Ribasim.set_flow!(p.graph, NodeID(:FlowBoundary, 1), NodeID(:Basin, 2), 1.0)
Ribasim.set_objective_priority!(allocation_model, p, u, t, priority_idx)
Ribasim.set_initial_values!(allocation_model, p, u, t)

println(p.allocation.allocation_models[1].problem)

Min F_abs_user_demand[UserDemand #13] + F_abs_user_demand[UserDemand #6] + F_abs_user_demand[UserDemand #3]
Subject to
 source[(FlowBoundary #1, Basin #2)] : F[(FlowBoundary #1, Basin #2)] ≤ 0
 source_user[UserDemand #13] : F[(UserDemand #13, Terminal #10)] ≤ 0
 source_user[UserDemand #6] : F[(UserDemand #6, Basin #5)] ≤ 0
 source_user[UserDemand #3] : F[(UserDemand #3, Basin #2)] ≤ 0
 fractional_flow[(TabulatedRatingCurve #7, FractionalFlow #8)] : F[(TabulatedRatingCurve #7, FractionalFlow #8)] - 0.6 F[(Basin #5, TabulatedRatingCurve #7)] ≤ 0
 fractional_flow[(TabulatedRatingCurve #7, FractionalFlow #9)] : F[(TabulatedRatingCurve #7, FractionalFlow #9)] - 0.4 F[(Basin #5, TabulatedRatingCurve #7)] ≤ 0
 flow_buffer_outflow[TabulatedRatingCurve #7] : F_flow_buffer_out[TabulatedRatingCurve #7] ≤ 0
 flow_conservation[LinearResistance #4] : -F[(LinearResistance #4, Basin #5)] + F[(Basin #5, LinearResistance #4)] + F[(Basin #2, LinearResistance #4)] - F[(LinearResistance #4, Basin #2)] = 0
 flow_conservation[Basin #5] : -F[(Basin #5, UserDemand #6)] + F[(LinearResistance #4, Basin #5)] - F[(Basin #5, LinearResistance #4)] + F[(UserDemand #6, Basin #5)] - F[(Basin #5, TabulatedRatingCurve #7)] = 0
 flow_conservation[DiscreteControl #11] : 0 = 0
 flow_conservation[TabulatedRatingCurve #7] : -F[(TabulatedRatingCurve #7, FractionalFlow #8)] - F[(TabulatedRatingCurve #7, FractionalFlow #9)] + F[(Basin #5, TabulatedRatingCurve #7)] - F_flow_buffer_in[TabulatedRatingCurve #7] + F_flow_buffer_out[TabulatedRatingCurve #7] = 0
 flow_conservation[Basin #12] : F[(FractionalFlow #9, Basin #12)] - F[(Basin #12, UserDemand #13)] = 0
 flow_conservation[Terminal #10] : F[(UserDemand #13, Terminal #10)] + F[(FractionalFlow #8, Terminal #10)] = 0
 flow_conservation[FractionalFlow #9] : -F[(FractionalFlow #9, Basin #12)] + F[(TabulatedRatingCurve #7, FractionalFlow #9)] = 0
 flow_conservation[Basin #2] : F[(FlowBoundary #1, Basin #2)] - F[(Basin #2, UserDemand #3)] + F[(UserDemand #3, Basin #2)] - F[(Basin #2, LinearResistance #4)] + F[(LinearResistance #4, Basin #2)] = 0
 flow_conservation[FractionalFlow #8] : F[(TabulatedRatingCurve #7, FractionalFlow #8)] - F[(FractionalFlow #8, Terminal #10)] = 0
 F[(UserDemand #13, Terminal #10)] ≥ 0
 F[(FlowBoundary #1, Basin #2)] ≥ 0
 F[(Basin #2, UserDemand #3)] ≥ 0
 F[(UserDemand #3, Basin #2)] ≥ 0
 F[(FractionalFlow #9, Basin #12)] ≥ 0
 F[(TabulatedRatingCurve #7, FractionalFlow #8)] ≥ 0
 F[(Basin #5, UserDemand #6)] ≥ 0
 F[(FractionalFlow #8, Terminal #10)] ≥ 0
 F[(TabulatedRatingCurve #7, FractionalFlow #9)] ≥ 0
 F[(Basin #12, UserDemand #13)] ≥ 0
 F[(LinearResistance #4, Basin #5)] ≥ 0
 F[(Basin #5, LinearResistance #4)] ≥ 0
 F[(Basin #2, LinearResistance #4)] ≥ 0
 F[(LinearResistance #4, Basin #2)] ≥ 0
 F[(UserDemand #6, Basin #5)] ≥ 0
 F[(Basin #5, TabulatedRatingCurve #7)] ≥ 0
 F_flow_buffer_in[TabulatedRatingCurve #7] ≥ 0
 F_flow_buffer_out[TabulatedRatingCurve #7] ≥ 0
 abs_positive_user_demand[UserDemand #13] : -F[(Basin #12, UserDemand #13)] + F_abs_user_demand[UserDemand #13] ≥ 0
 abs_positive_user_demand[UserDemand #6] : -F[(Basin #5, UserDemand #6)] + F_abs_user_demand[UserDemand #6] ≥ 0
 abs_positive_user_demand[UserDemand #3] : -F[(Basin #2, UserDemand #3)] + F_abs_user_demand[UserDemand #3] ≥ 0
 abs_negative_user_demand[UserDemand #13] : F[(Basin #12, UserDemand #13)] + F_abs_user_demand[UserDemand #13] ≥ 0
 abs_negative_user_demand[UserDemand #6] : F[(Basin #5, UserDemand #6)] + F_abs_user_demand[UserDemand #6] ≥ 0
 abs_negative_user_demand[UserDemand #3] : F[(Basin #2, UserDemand #3)] + F_abs_user_demand[UserDemand #3] ≥ 0