Allocation

1 Introduction

In water resources management, allocation deals with distributing limited water resources among competing demands in a fair and efficient manner. This is particularly important in situations where water demands exceed available supply, requiring a systematic approach to prioritize and optimize water distribution across the entire network. In Ribasim we formulate this with linear programming using JuMP and solve the problem using the the HiGHS solver.

1.1 Allocation Networks

Allocation is the process of assigning an optimized allocated flow rate to demand nodes in the physical layer of the model based on information about sources and their source priorities, the different demand nodes over various demand priorities, constraints introduced by nodes, local water availability and graph topology. It solves a linearized version of the problem formulated for the physical layer a period of time ahead of the physical layer, while simultaneously optimizing allocated flows and steering certain structures to get the water from (preferred) source to demand. The physics is implemented in an implicit manner, rather like the implicit Euler numerical method.

The allocation problem is solved per subnetwork of the Ribasim model. A subnetwork is defined by all nodes that have the same subnetwork ID in the Node table. Each subnetwork is used to formulate a linear 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.

When we build a Ribasim model in which we want to utilize allocation, it can be computationally efficient and stable to divide a large allocation network into a primary network which connects to multiple subnetworks. The primary network should be the main water system and treats each connected subnetwork as a single demand node.

We can then first solve an allocation optimization for the primary network, that calculates how to optimally divide the available water to each subnetwork. With this allocated amount for a subnetwork, we can optimize how to divide this water within the subnetwork to the demand nodes.

Subnetworks are defined by grouping nodes that share the same subnetwork ID in the Node table. The primary network needs to be assigned with subnetwork ID 1. This will always represent the main water system. Below is an example of a primary network with multiple subnetworks.

Code 
import ribasim_testmodels
import matplotlib.pyplot as plt
from IPython.display import Markdown, display

for model_name, model_constructor in ribasim_testmodels.constructors.items():

    if model_name!= "medium_primary_secondary_network":
        continue

    if model_constructor.__doc__ is not None:
        display(Markdown(model_constructor.__doc__))

    model = model_constructor()
    fig, ax = plt.subplots(figsize=(6, 4))
    model.plot(ax)
    ax.axis("off")
    plt.show()
    plt.close(fig)

Secondary networks:

  • Can only connect to the primary network via Pump or Outlet nodes
  • Cannot connect directly to other secondary networks
  • The demand of a subnetwork to the primary network is the sum of all unmet demands within the subnetwork

For each subnetwork, a linear optimization problem is formulated.

1.2 Mathematical formulation

In Ribasim, we optimize and assign allocated flow rates to demand nodes based on:

  • Information about sources and their priorities
  • Different demand nodes and their priorities
  • Constraints from nodes, local water availability and network topology

The optimization problem aims to minimize the difference between requested demands and supplied amounts, subject to water balance constraints:

\[ \begin{aligned} \text{min}\quad & z = \sum_{i=0}^{N_d} E_i(Q, S) \\ \text{s.t.}\quad & \frac{dS_j}{dt} = \sum_{k=1}^{N_l} Q_{k}\\ \end{aligned} \]

where:

  • \(z\) is the total error
  • \(N_d\) is the number of demands within a priority [-]
  • \(E_i\) is a deviation error between allocated and demand as function of flow and storage
  • \(S_j\) is the storage of node \(j\) [m³]
  • \(Q_k\) is the volumetric flow rate on link \(k\), which is positive for inflow and negative for outflow [m³/s]
  • \(N_l\) is the number of flow links connected to node \(j\)
  • \(t\) is the time [s]

1.3 Discretization

In the optimization problem, the water balance must be discretized over the allocation timestep.

\[\Delta t = t^{n+1} - t^{n}\]

Where the superscript \(^n\) denotes evaluation at the beginning of the timestep and \(^{n+1}\) at the end. We then use a backward Euler approximation:

\[\frac{dS}{dt} \approx \frac{S^{n+1} - S^{n}}{\Delta t} = \sum_{k=1}^{N_l} Q_{k}^{n+1} \]

Under the assumption that

  • flows are constant over \(\Delta t\)
  • coefficients in linearized equations are evaluated at the beginning of the timestep
  • The allocation timestep \(\Delta t\) is short enough that linearization remains valid

1.4 Linearization

To use efficient linear programming algorithms, nonlinear relationships can be linearized around the current time \(t^n\). For example, a nonlinear relationship between flow \(Q\), an upstream variable (\(x_1\)) and a downstream variable (\(x_2\)), can be linearized as:

\[Q(x_1,x_2) \approx Q(x_1^n,x_2^n) + \frac{\partial Q}{\partial x_1}(x_1^n,x_2^n)(x_1^{n+1}-x_1^n) + \frac{\partial Q}{\partial x_2}(x_1^n,x_2^n)(x_2^{n+1}-x_2^n)\]

This yields a linear programming problem that can be solved efficiently while maintaining acceptable accuracy for small changes around the operating point.

1.5 The high level algorithm

The allocation algorithm contains 2 types of optimization:

  • Demand collection, where water is allocated in the secondary networks with the sole purpose of finding out what the demand of the secondary network is from the primary network;
  • Allocation, where water is allocated in all subnetworks, and the amount of water that is allocated to demands is written as output and communicated to the physical layer.

The full algorithm goes through the following steps:

  1. Perform demand collection in the secondary subnetworks;
  2. Perform allocation in the primary network;
  3. Perform sequentially allocation (in the future maybe in parallel) in the subnetworks.

If no primary network is present, then step 1 and 2 are skipped and the entire network is treated as a single subnetwork.

2 Implementation

2.1 The basin profile

The level in a basin is generally a non-linear function of the storage (except when the area is constant). To incorporate this into a linear programming framework, we linearize the basin profile at the current timestep.

\[h^{n+1} \approx h^n + \frac{dh}{dS}(S^n)(S^{n+1} - S^n)\]

with

\[\frac{dh}{dS}(S^n) = \frac{1}{A^n}\]

and \(A^n\) is the area evaluated at the beginning of the time step

\[h^{n+1} \approx h^{n} + \frac{1}{A^n}(S^{n+1} - S^n)\]

2.2 Connector nodes

In Ribasim, connector nodes can determine the flow as function of upstream and/or downstream level. These are LinearResistance, ManningResistance and TabulatedRatingCurve nodes. Hence if we apply multi variable linearisation:

\[ Q^{n+1} = Q^{n} + \frac{\partial Q}{\partial h_1}(h_1^n,h_2^n)(h_1^{n+1} - h_1^{n}) + \frac{\partial Q}{\partial h_2}(h_1^n,h_2^n)(h_2^{n+1} - h_2^{n}) \]

However, we need to relate the level to a storage (in case of a basin). So we can substitute the linearised basin profile. For example if a non linear flow node connects downstream to a basin and upstream to a level boundary:

\[ Q_m^{n+1} = Q_m^{n} + \frac{1}{A^n}\frac{\partial Q_m}{\partial h_1}(h_1^n, h_2^n) (S_1^{n+1} - S_1^n) + \frac{\partial Q}{\partial h_2}(h_1^n,h_2^n)(h_2^{n+1} - h_2^{n}) \]

2.3 Basin forcings

2.4 Boundary nodes

We have the following boundary nodes in Ribasim:

  • The Terminal, where water simply leaves the model
  • The LevelBoundary, which yields fixed water levels; \(h^{n+1}_{lb}\) is read from the interpolated timeseries \(h_{lb}\) of the node or is a constant value
  • The FlowBoundary, which specifies a flow rate. Here the average outflow of the FlowBoundary in the physical layer over the previous \(\Delta t\) is used as a prediction of the flow rate in the next \(\Delta t\) over which the optimization takes place.
Note

The flow rate of a FlowBoundary is given as a timeseries so we could use the interpolation of that timeseries to compute the average flow in the coming \(\Delta t_\text{allocation}\). However, that would not always be accurate, since FlowBoundary nodes can be deactivated by Discrete Control which can not (easily) be anticipated.

3 Objectives (goals)

In Ribasim we have LinearResistance and ManningResistance nodes. The linear resistance node is already linear and therefore can be directly incorporated in the allocation problem. For ManningResistance we linearize the flow relationship with respect to the upstream and downstream level at the beginning of the allocation timestep.

For the resistance nodes the phantom storage introduced above is very helpful, because it prevents that according to the formula water has to flow from an empty basin with a high bottom to a lower basin.

3.1 Objectives (goals)

The allocation algorithm optimizes for a sequence of objectives in succession, adding constraints after optimization for each goal to retain the result of that optimization in the subsequent optimizations. This approach is known as goal programming.

3.1.1 Demand objectives

There are several types of demand nodes:

  • UserDemand (\(UD \subset S\)), which can have inflow demands for various demand priorities, and can consume a fraction of its abstraction and releases the rest back to the model;
  • LevelDemand (\(LD \subset S\)), which can have several level demands for various demand priorities, where higher levels must have lower priorities
  • FlowDemand (\(FD \subset S\)), which can have several inflow demands for various demand priorities, where no flow is consumed. The FlowDemand node gives this demand to another node

Demand objectives come in 2 categories:

  • flow demand (UserDemand, FlowDemand)
  • Level demand (LevelDemand)

We separate these out because these different types optimize for different quantities (flow and storage respectively), and combining demands of these two types within the same subnetwork within the same demand priority is not allowed. Given a demand priority \(d\), the objective function looks like

3.1.2 Formulation

In the following we assume lexicographic multi-objective optimization.

3.1.2.1 UserDemand and FlowDemand

Say we have a subnetwork with UserDemand nodes \(UD\) and FlowDemand nodes \(FD\) and each node \(v \in UD \cup FD\) has demands for priorities \(P_v\). Then for every one of these demand nodes for every priority for which they have a demand (at some point in time) we define a special flow decision variable which denotes how much is allocated to that node for that priority:

\[ 0 \le F^p_v \le d_v^p, \quad \forall v \in (UD \cup FD),\\; p \in P_v. \]

As indicated, each of these allocated flow variables is bounded between 0 and the corresponding demand (NOTE: in this text we assume that the demands are already expressed in the scaled flow unit). For UserDemand nodes we have that the sum of these allocated flow variables over the priorities is equal to the total flow into that node:

\[ F_{(b_\text{upstream}, v)} = \sum_{p \in P_v} F^p_v, \quad \forall v \in UD. \]

Note that this means that we enforce that only allocated flow can enter a UserDemand node, and so there can not be more flow into the UserDemand node than the total demand of that node. For FlowDemand we do allow that there is more flow through the node with the flow demand than the total demand (or even negative flow). To allow this freedom, we introduce a special ‘demand priority 0’ variable

\[ 0 \le F^0_v \le \frac{\text{MAX\_ABS\_FLOW}}{\text{scaling.flow}}, \quad \forall v \in FD \]

We only let \(F^0_v\) account for surplus flow, because if we allow that variable to be negative it can be used to allocate to demands. To make sure we still account for the possibility that the flow through the node with the flow demand is negative, we allow the allocated flow variable with the earliest priority to be negative, and account for this when writing the results.

With this we formulate the constraints

\[ F_{(b_\text{upstream}, v)} = \sum_{p} F^p_v, \quad \forall v \in FD. \]

We define relative error variables per demand node per priority for which they have a demand:

\[ E^p_v \ge 0, \quad \forall v \in (UD \cup FD),\\; p \in P_v. \]

The constraints on these relative error terms are similar as before:

\[ d_v^p \cdot E_v^p \ge d_v^p - F_v^p, \quad \forall v \in (UD \cup FD),\\; p \in P_v. \]

With these relative error terms we formulate an objective per demand priority for which any of these nodes have a demand (\(P^\text{flow} = \bigcup_{v \in UD \cup FD} P_v\)):

\[ \min \sum_{v,p} d_v^p \cdot E_v^p, \quad \forall p \in P^\text{flow}. \]

Note the multiplication by the demand here, which makes this effectively a minimization in the absolute error. There is a reason we formulate relative in stead of absolute errors: to optimize for a fair distribution of water. To that end we can formulate per demand priority a global relative allocation error:

\[ G^p = \frac{\sum_{v,p} d_v^p \cdot E_v^p}{\sum_{v,p} d^p_v}, \quad \forall p \in P^\text{flow}, \]

where in practice we multiply both sides by the denominator to prevent division by zero problems. Also note that the enumerator is precisely the expression for the aforementioned objective. Now in a subsequent objective, we want to penalize relative errors being larger than \(G^p\). So we define new errors

\[ \overline{E}_v^p \ge 0, \quad \forall v \in (UD \cup FD),\\; p \in P_v, \]

on which we define the constraints

\[ \overline{E}_v^p \ge E_v^p - G^p, \quad \forall v \in (UD \cup FD),\\; p \in P_v \]

and with which we formulate the objectives

\[ \min \sum_{v,p} \overline{E}_v^p, \quad \forall p \in P^\text{flow}. \]

Of course the order of objectives is important. For each demand priority of the flow type described here, first the weighted sum of the errors \(E_v^p\) must be minimized, then the sum of the errors \(\overline{E}_v^p\).

3.1.2.2 LevelDemand

Say we have Basins with LevelDemand \(B\) where each Basin \(b \in B\) has a minimum level \(h_{b,\min}^p\) for demand priorities \(P^\min_b\) and a maximum level \(h_{b,\max}^p\) for demand priorities \(P_b^\max\). Furthermore each Basin has a level to storage function \(s_b(h)\). For all Basins \(b\) in the subnetwork (not just the ones with a level demand) we also have

  • The known initial level at the beginning of the allocation timestep \(h_b^\text{init}\)
  • The decision variable of the change in storage over the allocation timestep \(\Delta S_b\)

For Basins with a level demand \(b \in B\) we additionally have

  • The lower absolute storage error

\[ E_{b, \text{lower}}^p \ge 0 \quad \forall b \in B, \\; p \in P^\text{min}_b \]

  • The incoming storage demand

\[ d^p_{b, +} = s\left(h_{b, \min}^p\right) - s_b\left(\text{clamp}\left(h_b^\text{init}, h_{b,\min}^{\text{prev}(p)}, h_{b,\min}^p \right)\right) \quad \forall b \in B, \\; p \in P^\text{min}_b \]

which is the volume of water needed to get to the minimum level of the current demand priority from the minimum level of the previous demand priority (or the starting level if that’s higher). If there is no priority in \(P_b^\min\) lower than the current priority \(p\) we say \(h^{\text{prev}(p)}_{b, \min}\) is equal to the Basin bottom.

  • The upper absolute storage error

\[ E_{b, \text{upper}}^p \ge 0 \quad \forall b \in B, \\; p \in P_b^\max \]

  • The outgoing storage demand:

\[ d^p_{b, -} = s_b\left(\text{clamp}\left(h_b^\text{init}, h_{b, \max}^p, h_{b, \max}^{\text{prev}(p)}\right)\right) - s\left(h_{b, \max}^p\right) \quad \forall b \in B, \\; p \in P^\text{max}_b \]

which is the volume of water needed to get from the maximum level of the current demand priority to the maximum level of the previous demand priority (or the starting level if that’s lower). If there is no priority in \(P_b^\max\) lower than the current priority \(p\) we say \(h^{\text{prev}(p)}_{b, \max} := \infty\).

We define the following constraints on the errors:

\[ E^p_{b, \text{lower}} \ge s\left(h^p_{b, \min}\right) -\left(s\left(h_b^\text{init}\right) + \Delta S_b\right) \quad \forall b \in S, \\; p \in P^\min_b \]

\[ E^p_{b, \text{upper}} \ge \left(s\left(h_b^\text{init}\right) + \Delta S_b\right) - s\left(h^p_{b, \max}\right) \quad \forall b \in S, \\; p \in P^\max_b. \]

Note that unlike with UserDemand and LevelDemand, the allocated amounts are not present as variables in the problem. The allocated amounts have to be derived from \(\Delta S_b\). Also note that the demands, which depend on the state of the physical layer and the start of the allocation timestep, do not appear in the LP problem. These are only computed for output.

Then given all this we define the objective for minimizing the absolute error sum per priority for which there is at least one level demand (\(P_\text{lower}^\text{level} = \bigcup_{b \in B} P^\min_b\), \(P_\text{upper}^\text{level} = \bigcup_{b \in B} P^\max_b\), \(P^\text{level} = P_\text{lower}^\text{level} \cup P_\text{upper}^\text{level}\)):

\[ \min \sum_{b; p\in P^\min_b} E^p_{b, \text{lower}} + \sum_{b \in B \\; : \\; p\in P^\max_b} E^p_{b, \text{upper}}, \qquad \forall p \in P^\text{level} \]

Now, just as in the previous section, we want to add a subsequent objective for fair distribution. Given the Basin area \(A_b\) (which is assumed to be constant over the allocation timestep), we can define the average lower and upper level error

\[ \Delta H^p_\text{lower} = \frac{\sum_{b \in B \\; : \\; p\in P^\min_b} E^p_{b, \text{lower}}}{\sum_{b \in B \\; : \\; p\in P^\min_b} A_b}, \qquad \forall p \in P_\text{lower}^\text{level}; \]

\[ \Delta H^p_\text{upper} = \frac{\sum_{b \in B \\; : \\; p\in P^\max_b} E^p_{b, \text{upper}}}{\sum_{b \in B \\; : \\; p\in P^\max_b} A_b}, \qquad \forall p \in P_\text{upper}^\text{level}; \]

Again we define new errors:

\[ \overline{E}_{b, \text{lower}}^p \ge 0, \quad \forall b \in B, \\; p \in P_b^\min; \]

\[ \overline{E}_{b, \text{upper}}^p \ge 0, \quad \forall b \in B, \\; p \in P_b^\max. \]

On these errors we define the following constraints:

\[ \overline{E}^p_{b, \text{lower}} \ge \frac{E^p_{b, \text{lower}}}{A_b} - \Delta H^p_\text{lower}, \quad \forall b \in B, \\;p \in P_b^\text{min}; \]

\[ \overline{E}^p_{b, \text{upper}} \ge \frac{E^p_{b, \text{upper}}}{A_b} - \Delta H^p_\text{upper}, \quad \forall b \in B, \\;p \in P_b^\text{max}, \]

and the following objectives:

\[ \min \sum_{b \in B \\; : \\; p \in P_b^\min} \overline{E}^p_{b, \text{lower}} + \sum_{b \in B \\; : \\; p \in P_b^\max} \overline{E}^p_{b, \text{upper}}, \qquad \forall p \in P^\text{level}. \]

3.2 Control

3.2.1 Control by allocation

When Pumps and Outlets are part of a subnetwork, they can be controlled by allocation. To accomplish this, they must be given the special control state Ribasim.allocation in the static table for Pump or Outlet. When a Pump or Outlet has this control state, the flow through the node is only bounded by its capacity and Pump or Outlet specific constraints based on the difference between upstream and downstream level (see Pump equations or Outlet equations). After all goals have been optimized for, the flow rate through the Pump or Outlet is communicated to the physical layer.

3.2.2 Interaction of allocation with other control systems

There are several other control systems in Ribasim:

  • DiscreteControl: This node type can change parameters of other nodes in between time steps of the solver of the physical layer. If the affected node is within a subnetwork, the parameter change will also be taken into account in the next allocation run. So the parameters in the allocation layer are always up to date with those in the physical layer, but the allocation algorithm cannot anticipate parameter changes from DiscreteControl that occur within the allocation time step which is being optimized over, so this can be a source of discrepancies between the physical layer and the allocation layer.

  • ContinuousControl: The continuous nature of this control type is not taken into account. The flow rate the controlled structure has at the start of the allocation optimization step will be extrapolated in a constant manner.

  • PidControl: The continuous nature of this control type is not taken into account. The flow rate the controlled structure has at the start of the allocation optimization step will be extrapolated in a constant manner.

The following data of the parameters and state of a Ribasim model are relevant for the allocation problem.

3.3 Schematisation input

3.3.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 each other. Nodes can also not be part of any subnetwork.

3.3.2 Source flows

Sources are indicated by a set of links in the subnetwork \[ E_S^\text{source} \subset E, \] which are automatically inferred as all links 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 flow over this link \[ \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 links are either coming from a boundary/source node (e.g. a level or flow boundary) or connect the primary network to a subnetwork. For the definition of \(Q_{ij}\) see the formal model description.

3.3.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.3.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 through node \(i\) at all priorities \(p < p_{\text{fd}}\).

3.4 Simulation (physical layer) input

3.4.1 Constraining factors

3.4.1.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.4.1.2 UserDemand return flows

UserDemand nodes dictate proportional relationships between flows over links in the subnetwork. The return factor is given by \(0 \le r_i(t) \le 1, i \in U_S\).

3.5 The subnetwork

The subnetwork consists of a set of nodes \(S \subset V\) and links

\[ E_S = (S \times S) \cup E_S^\text{source}, \]

i.e. the links that lie within the subnetwork together with the source links (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.5.1 Capacities

Each link 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 link capacity is infinite if there is nothing in the model constraining the capacity.

The capacities are determined in different ways:

  • If an link 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 link \(e \in E_S\) is given by the smallest max_flow_rate of the nodes along the equivalent links in the subnetwork. If there are no nodes with a max_flow_rate, the link capacity is infinite. If the max_flow_rate is time-dependent, only the value at the starttime of the simulation is considered;
  • If the link is a source, the capacity of the link is given by the flow rate of that source;

There are also capacities for special links:

  • \(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.