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.

#TODO: Add model image with primary network and secondary networks.

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:

Perform demand collection in the secondary subnetworks;
Perform allocation in the primary network;
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.

At the time of writing it is not completely clear yet which parts of the physics will be added as hard constraints and which ones will be optimized for (and how).

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

#TODO: Write out objective functions.

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. When a Pump or Outlet has this control state, the flow trough the node is only bounded by its capacity and #TODO: Pump or Node specific constraints based on the difference between upstream and downstream level.. After all goals have been optimized for, the flow rate trough 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 #TODO: This has not been implemented yet.. 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: #TODO: Undecided. Some ContinuousControl relationships could be implemented in allocation, but this requires validation.
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.

See here.

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 trough 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.