Avoid unconstrained sections in allocation networks (#1429)

Fixes #1208.

I fixed the `UserDemand` point more generally: no sources other than the
main network should be available when collecting subnetwork demands. We
do allow multiple inlets from the main network to a subnetwork. We do
not test this yet, but this also allows for insufficiently constraint
distribution of the subnetwork demands over these inlets. A preference
ordering over sources (#565)
could resolve this.

core/src/allocation_init.jl

@@ -164,8 +164,9 @@ function add_variables_basin!(
 end
 
 """
-Add the variables for supply/demand of the buffer of a node with a flow demand to the problem.
-The variable indices are the node IDs of the nodes with a flow demand in the subnetwork.
+Add the variables for supply/demand of the buffer of a node with a flow demand
+or fractional flow outneighbors to the problem.
+The variable indices are the node IDs of the nodes with a buffer in the subnetwork.
 """
 function add_variables_flow_buffer!(
     problem::JuMP.Model,
@@ -175,9 +176,11 @@ function add_variables_flow_buffer!(
     (; graph) = p
 
     # Collect the nodes in the subnetwork that have a flow demand
+    # or fractional flow outneighbors
     node_ids_flow_demand = NodeID[]
     for node_id in graph[].node_ids[subnetwork_id]
-        if has_external_demand(graph, node_id, :flow_demand)[1]
+        if has_external_demand(graph, node_id, :flow_demand)[1] ||
+           has_fractional_flow_outneighbors(graph, node_id)
             push!(node_ids_flow_demand, node_id)
         end
     end
@@ -366,13 +369,7 @@ function add_constraints_conservation_node!(
         is_source_sink = node_id.type in
         [NodeType.FlowBoundary, NodeType.LevelBoundary, NodeType.UserDemand]
 
-        # No flow conservation on nodes with FractionalFlow outneighbors
-        has_fractional_flow_outneighbors = any(
-            outflow_id.type == NodeType.FractionalFlow for
-            outflow_id in outflow_ids(graph, node_id)
-        )
-
-        if is_source_sink | has_fractional_flow_outneighbors
+        if is_source_sink
             continue
         end
 
@@ -399,8 +396,9 @@ function add_constraints_conservation_node!(
             push!(outflows_node, F_basin_in[node_id])
         end
 
-        # If the node has a flow demand
-        if has_external_demand(graph, node_id, :flow_demand)[1]
+        # If the node has a buffer
+        if has_external_demand(graph, node_id, :flow_demand)[1] ||
+           has_fractional_flow_outneighbors(graph, node_id)
             push!(inflows_node, F_flow_buffer_out[node_id])
             push!(outflows_node, F_flow_buffer_in[node_id])
         end

core/src/allocation_optim.jl

@@ -247,7 +247,7 @@ end
 
 """
 Set the capacities of the sources in the subnetwork
-as the latest instantaneous flow out of the source in the physical layer
+as the average flow over the last Δt_allocation of the source in the physical layer
 """
 function set_initial_capacities_source!(
     allocation_model::AllocationModel,
@@ -996,6 +996,26 @@ function set_initial_values!(
     return nothing
 end
 
+"""
+Set the capacities of all edges that denote a source to 0.0.
+"""
+function empty_sources!(allocation_model::AllocationModel, allocation::Allocation)::Nothing
+    (; problem) = allocation_model
+    (; subnetwork_demands) = allocation
+
+    for constraint_set_name in [:source, :source_user, :basin_outflow, :flow_buffer_outflow]
+        constraint_set = problem[constraint_set_name]
+        for key in only(constraint_set.axes)
+            # Do not set the capacity to 0.0 if the edge
+            # is a main to subnetwork connection edge
+            if key ∉ keys(subnetwork_demands)
+                JuMP.set_normalized_rhs(constraint_set[key], 0.0)
+            end
+        end
+    end
+    return nothing
+end
+
 """
 Update the allocation optimization problem for the given subnetwork with the problem state
 and flows, solve the allocation problem and assign the results to the UserDemand.
@@ -1027,7 +1047,11 @@ function allocate!(
 
     set_initial_capacities_inlet!(allocation_model, p, optimization_type)
 
-    if optimization_type != OptimizationType.collect_demands
+    if optimization_type == OptimizationType.collect_demands
+        # When collecting demands, only flow should be available
+        # from the main to subnetwork connections
+        empty_sources!(allocation_model, allocation)
+    else
         set_initial_values!(allocation_model, p, u, t)
     end
 

core/src/util.jl

@@ -731,3 +731,8 @@ function get_discrete_control_indices(discrete_control::DiscreteControl, conditi
         condition_idx_now += l
     end
 end
+
+has_fractional_flow_outneighbors(graph::MetaGraph, node_id::NodeID)::Bool = any(
+    outneighbor_id.type == NodeType.FractionalFlow for
+    outneighbor_id in outflow_ids(graph, node_id)
+)

core/test/allocation_test.jl

@@ -271,7 +271,7 @@ end
     @test user_demand.allocated[7, :] ≈ [0.001, 0.0, 0.0]
 end
 
-@testitem "subnetworks with sources" begin
+@testitem "Subnetworks with sources" begin
     using SQLite
     using Ribasim: NodeID, OptimizationType
     using ComponentArrays: ComponentVector
@@ -307,7 +307,7 @@ end
     # See the difference between these values here and in
     # "allocation with main network optimization problem", internal sources
     # lower the subnetwork demands
-    @test subnetwork_demands[(NodeID(:Basin, 2), NodeID(:Pump, 11))] ≈ [3.1, 4.0, 0.0]
+    @test subnetwork_demands[(NodeID(:Basin, 2), NodeID(:Pump, 11))] ≈ [4.0, 4.0, 0.0]
     @test subnetwork_demands[(NodeID(:Basin, 6), NodeID(:Pump, 24))] ≈ [0.004, 0.0, 0.0]
     @test subnetwork_demands[(NodeID(:Basin, 10), NodeID(:Pump, 38))][1:2] ≈ [0.001, 0.001]
 end

docs/core/allocation.qmd

@@ -42,7 +42,7 @@ 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 allocation interval $\Delta t_{\text{alloc}}$ of the of the flow over this edge
+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'
 $$
@@ -75,7 +75,7 @@ for all $i \in FD_S$. Here $d^{p_{\text{df}}}$ is given by the original flow dem
 
 ### 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 allocation interval $\Delta t_{\text{alloc}}$ of the vertical fluxes (precipitation, evaporation, infiltration and drainage) for each basin:
+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.
 $$
@@ -138,8 +138,8 @@ The optimization problem for a subnetwork is a linear optimization problem consi
 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 respectively;
-- The flows $F^\text{buffer out}_{i}, F^\text{buffer in}_{i} \ge 0$ for all $i \in FD_S$ supplied and consumed by the flow buffers of nodes with a flow demand respectively.
+- 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.
 
 ## The optimization objective
 
@@ -183,19 +183,18 @@ For convenience, we use the notation
 
 for the set of in-neighbors and out-neighbors of a node in the network respectively.
 
-- Flow conservation: For the basins in the allocation network we have that
+- 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:
 $$
-    F^\text{basin in}_k + \sum_{j \in V^{\text{out}}_S(k)} F_{kj} = F^\text{basin out}_k + \sum_{i \in V^{\text{in}}_S(k)} F_{ik}, \quad \forall k \in B_S .
+    \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.
 $$  {#eq-flowconservationconstraintbasin}
-We have the same constraint without the basin terms for nodes that have flow edges as inneighbors (except if this node also happens to be a basin).
-For nodes which have a flow demand we have
-$$
-F_{kj} + F^\text{buffer in}_k  = F^\text{flow in}_k + F_{ik}, \quad \forall k \in FD_S, \quad V^{\text{in}}_S(k) = \{i\},\;
-     V^{\text{out}}_S(k) = \{j\}.
-$$ {#eq-flowconservationconstraintflowdemand}
+
+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](#the-optimization-variables)) there is a special outflow $F^{\text{buffer in}}_k$ and a special inflow $F^{\text{buffer out}}_k$.
 
 :::{.callout-note}
-In @eq-flowconservationconstraintbasin and @eq-flowconservationconstraintflowdemand, 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.
+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:
@@ -209,7 +208,7 @@ $$
 $$
 
 :::{.callout-note}
-When performing subnetwork demand collection, these capacities are set to $\infty$ for edges which connect the main network to a subnetwork.
+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: