Allocation

1 Introduction

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.

There can be a special subnetwork called primary network (with subnetwork ID 1), which represents the main or primary water system. Other subnetworks are called secondary networks and are only allowed to connect to the primary network via a Pump or Outlet; secondary networks are not allowed to be connected to each other. Secondary networks automatically formulate demands to the primary network if there is a shortage of water within the subnetwork to fulfill all demands.

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

2 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 allocation for the other subnetworks.

If no primary network is present, then step 1 and 2 are skipped.

3 The optimization problems for the subnetworks

3.1 Water balance

Just as the physical layer, the allocation layer must satisfy the water balance, in this case over the allocation timestep \(\Delta t_\text{allocation}\). To explain how this is done, we give an algebraic description of some basic aspects of the subnetworks. Say we denote the set of all nodes in the model by \(V\), and the subset that forms a subnetwork by \(S \subseteq V\). We also have the the set of all flow links in the model \(L \subseteq V \times V\), i.e. the links are ordered pairs of nodes. The set of links associated with the subnetwork is derived from \(S\) as

\[ E_S = \{l \in L \; : \; l \cap S \ne\emptyset\}, \]

#TODO: Add image with which links are assumed to be part of the subnetwork.

i.e. all links in the model that have at least one node in the subnetwork. For each link in the subnetwork we define a flow variable (\(\text{m}^3/\text{s}\)):

\[ -F_\max \le \mathcal{F}_l \le F_\max, \quad l \in E_S. \]

Here each \(\mathcal{F}_l\) denotes the flow over the link \(l\) over the timestep. These are decision variables, i.e. values that are an output of the optimization, which we denote in cursive script. The variable can attain values between the bounds \(-F_\max\) and \(F_\max\), where \(F_\max = 500.000 \;\text{m}^3/\text{s}\). The bound \(F_\max\) is chosen as a reasonable maximum value that flows can attain. We assume these flows to be constant over the period \(\Delta t_\text{allocation}\). We will see later that certain parameters of the model can put stricter bounds on these flows.

For the basins \(B_S \subset S\) in the subnetwork we define storage parameters (\(\text{m}^3\)) for the start of the allocation time step \(\Delta t_\text{allocation}\):

\[ S^\text{start}_b, \quad b \in B_S. \]

These are set to the values that hey have in the physical layer at the start of the allocation timestep. The storages at the end of this timestep are decision variables:

\[ 0 \le \mathcal{S}^\text{end}_b \le S_b^\max, b \in B_S, \]

#TODO: We might decide to relax this upper bound a bit.

where \(S_b^\max\) is the largest storage within the Basin profile. Lastly we have a horizontal forcing variable for each basin

\[ \mathcal{H}^\text{end}, \quad b \in B_S, \]

whose value will be discussed later. Given these variables, we can set up a water balance constraint for each basin in the subnetwork:

\[ \mathcal{S}_b^\text{end} - S_b^\text{start} = \Delta t_\text{allocation} \left( \mathcal{H}_b + \sum_{(n, b) \in E_S} \mathcal{F}_{(n,b)} - \sum_{(b, n) \in E_S} \mathcal{F}_{(b,n)} \right), \quad b \in B_S. \]

3.2 The Basin profile

For most of the physical layer, flows depend on levels. Therefore level variables have to be introduced:

\[ L_b^\min \ge \mathcal{L}^\text{end}_b \ge L_b^\max, \quad b \in B_S, \]

#TODO: We might decide to relax this upper bound a bit, in accordance with the relaxation of the upper bound of the storage.

where \(L_b^\min\) is the bottom level of the Basin and \(L_b^\max\) is the highest level available in the Basin profile. Note that this means that the problem becomes infeasible if the level in a basin exceeds this value, and the solver will attempt to avoid this from happening. The level in a basin depends on the storage via the Basin profile, here denoted by the function \(P_b\):

\[ \mathcal{L}^\text{end}_b = P_b(\mathcal{S}^\text{end}_b) \]

#TODO: Explain how the level(storage) relationship is linearized, and how some ‘phantom storage’ is introduced into the profile.

Just like the initial storages, each Basin has an initial level \(L_b^\text{start}\) known from the physical layer which we need to compute flows.

3.3 Physical processes

3.3.1 Basin forcings

3.3.2 Boundary nodes

We have the following boundary nodes in Ribasim:

The Terminal, where water can leave the model;
The LevelBoundary, which yields fixed water levels; \(L^\text{end}_\text{fb} = h_{fb}(t + \Delta t_\text{allocation})\) given the interpolated timeseries \(h_{fb}\) of the node (or 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_\text{allocation}\) is used as a prediction of the flow rate in the next \(\Delta t_\text{allocation}\) 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.3.3 Tabulated Rating Curve

Since the \(Q(h)\) relationships for TabulatedRatingCurve nodes (see here) are already linear, these can be directly incorporated into the allocation problem.

3.3.4 Resistance nodes

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.4 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.4.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.5 Control

3.5.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.5.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.

4 Output

See here.