Add technical documentation for allocation (#1579)

Fixes #1540 We now have a technical documentation for allocation! --------- Co-authored-by: Bart de Koning <[email protected]> Co-authored-by: Bart de Koning <[email protected]>
Deltares · Jul 31, 2024 · 38cab4f · 38cab4f
1 parent 900afec
commit 38cab4f
Show file tree

Hide file tree

Showing 6 changed files with 380 additions and 26 deletions.
diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl
@@ -3,7 +3,7 @@ function find_subnetwork_connections!(p::Parameters)::Nothing
     (; allocation, graph, allocation) = p
     n_priorities = length(allocation.priorities)
     (; subnetwork_demands, subnetwork_allocateds) = allocation
-    # Find edges where the source node has subnetwork id 1 and the
+    # Find edges (node_id, outflow_id) where the source node has subnetwork id 1 and the
     # destination node subnetwork id ≠1
     for node_id in graph[].node_ids[1]
         for outflow_id in outflow_ids(graph, node_id)
@@ -22,8 +22,25 @@ function find_subnetwork_connections!(p::Parameters)::Nothing
     return nothing
 end
 
+function get_main_network_connections(
+    p::Parameters,
+    subnetwork_id::Int32,
+)::Vector{Tuple{NodeID, NodeID}}
+    (; allocation) = p
+    (; subnetwork_ids, main_network_connections) = allocation
+    idx = findsorted(subnetwork_ids, subnetwork_id)
+    if isnothing(idx)
+        error("Invalid allocation network ID $subnetwork_id.")
+    else
+        return main_network_connections[idx]
+    end
+    return
+end
+
 """
-Get the fixed capacity of the edges in the subnetwork
+Get the fixed capacity (∈[0,∞]) of the edges in the subnetwork in a JuMP.Containers.SparseAxisArray,
+which is a type of sparse arrays that in this case takes NodeID in stead of Int as indices.
+E.g. capacity[(node_a, node_b)] gives the capacity of edge (node_a, node_b).
 """
 function get_subnetwork_capacity(
     p::Parameters,
@@ -61,6 +78,7 @@ function get_subnetwork_capacity(
                 capacity_edge = min(capacity_edge, capacity_node_dst)
             end
 
+            # Set the capacity
             capacity[edge_metadata.edge] = capacity_edge
 
             # If allowed by the nodes from this edge,
@@ -82,7 +100,7 @@ const allocation_source_nodetypes =
 
 """
 Add the edges connecting the main network work to a subnetwork to both the main network
-and subnetwork allocation network.
+and subnetwork allocation network (defined by their capacity objects).
 """
 function add_subnetwork_connections!(
     capacity::JuMP.Containers.SparseAxisArray{Float64, 2, Tuple{NodeID, NodeID}},
@@ -151,6 +169,8 @@ function add_variables_basin!(
     subnetwork_id::Int32,
 )::Nothing
     (; graph) = p
+
+    # Get the node IDs from the subnetwork for basins that have a level demand
     node_ids_basin = [
         node_id for
         node_id in graph[].node_ids[subnetwork_id] if graph[node_id].type == :basin &&
@@ -313,7 +333,8 @@ function add_constraints_conservation_node!(
 
     for node_id in node_ids
 
-        # No flow conservation constraint on sources/sinks
+        # If a node is a source or a sink (i.e. a boundary node),
+        # there is no flow conservation on that node
         is_source_sink = node_id.type in
         [NodeType.FlowBoundary, NodeType.LevelBoundary, NodeType.UserDemand]
 
@@ -368,6 +389,56 @@ function add_constraints_conservation_node!(
     return nothing
 end
 
+"""
+Add the fractional flow constraints to the allocation problem.
+The constraint indices are allocation edges over a fractional flow node.
+
+Constraint:
+flow after fractional_flow node <= fraction * inflow
+"""
+function add_constraints_fractional_flow!(
+    problem::JuMP.Model,
+    p::Parameters,
+    subnetwork_id::Int32,
+)::Nothing
+    (; graph, fractional_flow) = p
+    F = problem[:F]
+    node_ids = graph[].node_ids[subnetwork_id]
+
+    # Find the nodes in this subnetwork with a FractionalFlow
+    # outneighbor, and collect the corresponding flow fractions
+    # and inflow variable
+    edges_to_fractional_flow = Tuple{NodeID, NodeID}[]
+    fractions = Dict{Tuple{NodeID, NodeID}, Float64}()
+    inflows = Dict{NodeID, JuMP.AffExpr}()
+
+    # Find edges of the form (node_id, outflow_id) where outflow_id
+    # is for a FractionalFlow node
+    for node_id in node_ids
+        for outflow_id in outflow_ids(graph, node_id)
+            if outflow_id.type == NodeType.FractionalFlow
+                edge = (node_id, outflow_id)
+                push!(edges_to_fractional_flow, edge)
+                fractions[edge] = fractional_flow.fraction[outflow_id.idx]
+                inflows[node_id] = sum([
+                    F[(inflow_id, node_id)] for inflow_id in inflow_ids(graph, node_id)
+                ])
+            end
+        end
+    end
+
+    # Create the constraints if there is at least one
+    if !isempty(edges_to_fractional_flow)
+        problem[:fractional_flow] = JuMP.@constraint(
+            problem,
+            [edge = edges_to_fractional_flow],
+            F[edge] <= fractions[edge] * inflows[edge[1]],
+            base_name = "fractional_flow"
+        )
+    end
+    return nothing
+end
+
 """
 Add the Basin flow constraints to the allocation problem.
 The constraint indices are the Basin node IDs.

diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl
@@ -40,10 +40,12 @@ function set_objective_priority!(
     (; main_network_connections, subnetwork_demands) = allocation
     F = problem[:F]
 
+    # Initialize an empty quadratic expression for the objective
     ex = JuMP.QuadExpr()
 
-    # Terms for subnetworks as UserDemand
+    # Terms for subnetworks acting as UserDemand on the main network
     if is_main_network(subnetwork_id)
+        # Loop over the connections between main and subnetwork
         for connections_subnetwork in main_network_connections[2:end]
             for connection in connections_subnetwork
                 d = subnetwork_demands[connection][priority_idx]
@@ -53,7 +55,7 @@ function set_objective_priority!(
         end
     end
 
-    # Terms for UserDemand nodes and LevelDemand nodes
+    # Terms for UserDemand nodes and FlowDemand nodes
     for edge in keys(capacity.data)
         to_node_id = edge[2]
 
@@ -91,6 +93,7 @@ function set_objective_priority!(
         add_objective_term!(ex, d, F_ld)
     end
 
+    # Add the new objective to the problem
     new_objective = JuMP.@expression(problem, ex)
     JuMP.@objective(problem, Min, new_objective)
     return nothing
@@ -608,6 +611,8 @@ function adjust_demands!(
     F = problem[:F]
 
     for node_id in flow_demand.node_id
+
+        # Only update data for FlowDemand nodes in the current subnetwork
         if graph[node_id].subnetwork_id != subnetwork_id
             continue
         end
@@ -653,6 +658,9 @@ function adjust_capacities_buffer!(allocation_model::AllocationModel)::Nothing
 
     for node_id in only(constraints_flow_buffer.axes)
         constraint = constraints_flow_buffer[node_id]
+
+        # The capacity should not be able to get below 0, but can get to small negative numbers,
+        # probably due to floating point errors. Therefore new_capacity = max(0, new_capacity) is applied.
         buffer_capacity = max(
             0.0,
             JuMP.normalized_rhs(constraint) + JuMP.value(F_flow_buffer_in[node_id]) -
@@ -667,6 +675,10 @@ end
 Set the capacity of the outflow edge from a node with a flow demand:
 - To Inf if the current priority is other than the priority of the flow demand
 - To 0.0 if the current priority is equal to the priority of the flow demand
+
+This is done so that flow can go towards the node with the flow demand into its buffer,
+to avoid the problem that the flow has nowhere to go after this node and to make sure
+that this flow can be used for later priorities
 """
 function set_capacities_flow_demand_outflow!(
     allocation_model::AllocationModel,
@@ -1074,10 +1086,6 @@ function allocate_demands!(
     (; subnetwork_id) = allocation_model
     (; priorities) = allocation
 
-    if is_main_network(subnetwork_id)
-        @assert optimization_type == OptimizationType.allocate "For the main network no demands have to be collected"
-    end
-
     set_initial_capacities_inlet!(allocation_model, p, optimization_type)
 
     set_initial_values!(allocation_model, p, u, t)

diff --git a/core/src/callback.jl b/core/src/callback.jl
@@ -366,21 +366,6 @@ function get_allocation_model(p::Parameters, subnetwork_id::Int32)::AllocationMo
     end
 end
 
-function get_main_network_connections(
-    p::Parameters,
-    subnetwork_id::Int32,
-)::Vector{Tuple{NodeID, NodeID}}
-    (; allocation) = p
-    (; subnetwork_ids, main_network_connections) = allocation
-    idx = findsorted(subnetwork_ids, subnetwork_id)
-    if isnothing(idx)
-        error("Invalid allocation network ID $subnetwork_id.")
-    else
-        return main_network_connections[idx]
-    end
-    return
-end
-
 function set_control_params!(p::Parameters, node_id::NodeID, control_state::String)::Nothing
     (; discrete_control, allocation) = p
     (; control_mappings) = discrete_control

diff --git a/core/src/validation.jl b/core/src/validation.jl
@@ -624,7 +624,9 @@ function valid_discrete_control(p::Parameters, config::Config)::Bool
 end
 
 """
-The source nodes must only have one allocation outneighbor and no allocation inneighbors.
+An allocation source edge is valid if either:
+    - The edge connects the main network to a subnetwork
+    - The edge comes from a source node
 """
 function valid_sources(
     p::Parameters,
@@ -635,12 +637,18 @@ function valid_sources(
 
     errors = false
 
+    # Loop over edges that were assigned a capacity
     for edge in keys(capacity.data)
+
+        # For an edge (id_a, id_b) in the physical model
+        # the reverse (id_b, id_a) can exist in the allocation subnetwork
         if !haskey(graph, edge...)
             edge = reverse(edge)
         end
 
         (id_source, id_dst) = edge
+
+        # Whether the current edge is a source for the current subnetwork
         if graph[edge...].subnetwork_id_source == subnetwork_id
             from_source_node = id_source.type in allocation_source_nodetypes
 

diff --git a/docs/_quarto.yml b/docs/_quarto.yml
@@ -63,6 +63,7 @@ website:
         - reference/index.qmd
         - reference/usage.qmd
         - reference/validation.qmd
+        - reference/allocation.qmd
         - reference/python/index.qmd
         - reference/test-models.qmd
         - section: "Nodes"
@@ -90,6 +91,7 @@ website:
         - dev/addnode.qmd
         - dev/python.qmd
         - dev/qgis.qmd
+        - dev/allocation.qmd
         - dev/bmi.qmd
         - dev/ci.qmd
         - dev/release.qmd