diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 6ef809bd8..6fd3ddbd5 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -295,30 +295,6 @@ end const allocation_source_nodetypes = Set{NodeType.T}([NodeType.LevelBoundary, NodeType.FlowBoundary]) -""" -Remove allocation UserDemand return flow edges that are upstream of the UserDemand itself. -""" -function avoid_using_own_returnflow!(p::Parameters, allocation_network_id::Int32)::Nothing - (; graph) = p - node_ids = graph[].node_ids[allocation_network_id] - edge_ids = graph[].edge_ids[allocation_network_id] - node_ids_user_demand = - [node_id for node_id in node_ids if node_id.type == NodeType.UserDemand] - - for node_id_user_demand in node_ids_user_demand - node_id_return_flow = only(outflow_ids_allocation(graph, node_id_user_demand)) - if allocation_path_exists_in_graph(graph, node_id_return_flow, node_id_user_demand) - edge_metadata = graph[node_id_user_demand, node_id_return_flow] - graph[node_id_user_demand, node_id_return_flow] = - @set edge_metadata.allocation_flow = false - empty!(edge_metadata.node_ids) - delete!(edge_ids, (node_id_user_demand, node_id_return_flow)) - @debug "The outflow of $node_id_user_demand is upstream of the UserDemand itself and thus ignored in allocation solves." - end - end - return nothing -end - """ Add the edges connecting the main network work to a subnetwork to both the main network and subnetwork allocation network. @@ -359,9 +335,6 @@ function allocation_graph( error("Errors in sources in allocation network.") end - # Discard UserDemand return flow in allocation if this leads to a closed loop of flow - avoid_using_own_returnflow!(p, allocation_network_id) - return capacity end @@ -512,6 +485,31 @@ function add_constraints_capacity!( return nothing end +""" +Add capacity constraints to the outflow edge of UserDemand nodes. +The constraint indices are the UserDemand node IDs. + +Constraint: +flow over UserDemand edge outflow edge <= cumulative return flow from previous priorities +""" +function add_constraints_user_source!( + problem::JuMP.Model, + p::Parameters, + allocation_network_id::Int32, +)::Nothing + (; graph) = p + F = problem[:F] + node_ids = graph[].node_ids[allocation_network_id] + node_ids_user = [node_id for node_id in node_ids if node_id.type == NodeType.UserDemand] + + problem[:source_user] = JuMP.@constraint( + problem, + [node_id = node_ids_user], + F[(node_id, outflow_id(graph, node_id))] <= 0.0 + ) + return nothing +end + """ Add the source constraints to the allocation problem. The actual threshold values will be set before each allocation solve. @@ -693,41 +691,6 @@ function add_constraints_conservation_basin!( return nothing end -""" -Add the UserDemand returnflow constraints to the allocation problem. -The constraint indices are UserDemand node IDs. -Constraint: -outflow from user_demand <= return factor * inflow to user_demand -""" -function add_constraints_user_demand_returnflow!( - problem::JuMP.Model, - p::Parameters, - allocation_network_id::Int32, -)::Nothing - (; graph, user_demand) = p - F = problem[:F] - - node_ids = graph[].node_ids[allocation_network_id] - node_ids_user_demand_with_returnflow = [ - node_id for node_id in node_ids if node_id.type == NodeType.UserDemand && - !isempty(outflow_ids_allocation(graph, node_id)) - ] - problem[:return_flow] = JuMP.@constraint( - problem, - [node_id_user_demand = node_ids_user_demand_with_returnflow], - F[( - node_id_user_demand, - only(outflow_ids_allocation(graph, node_id_user_demand)), - )] <= - user_demand.return_factor[findsorted(user_demand.node_id, node_id_user_demand)] * F[( - only(inflow_ids_allocation(graph, node_id_user_demand)), - node_id_user_demand, - )], - base_name = "return_flow", - ) - return nothing -end - """ Minimizing |expr| can be achieved by introducing a new variable expr_abs and posing the following constraints: @@ -971,7 +934,7 @@ function allocation_problem( add_constraints_capacity!(problem, capacity, p, allocation_network_id) add_constraints_source!(problem, p, allocation_network_id) - add_constraints_user_demand_returnflow!(problem, p, allocation_network_id) + add_constraints_user_source!(problem, p, allocation_network_id) add_constraints_fractional_flow!(problem, p, allocation_network_id) add_constraints_basin_flow!(problem) add_constraints_flow_demand_outflow!(problem, p, allocation_network_id) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 09fad5d3d..17e9eab37 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -424,6 +424,48 @@ function adjust_capacities_basin!( return nothing end +""" +Set the initial capacities of the UserDemand return flow sources to 0. +""" +function set_initial_capacities_returnflow!(allocation_model::AllocationModel)::Nothing + (; problem) = allocation_model + constraints_outflow = problem[:source_user] + + for node_id in only(constraints_outflow.axes) + constraint = constraints_outflow[node_id] + capacity = 0.0 + JuMP.set_normalized_rhs(constraint, capacity) + end + return nothing +end + +""" +Add the return flow fraction of the inflow to the UserDemand nodes +to the capacity of the outflow source. +""" +function adjust_capacities_returnflow!( + allocation_model::AllocationModel, + p::Parameters, +)::Nothing + (; graph, user_demand) = p + (; problem) = allocation_model + constraints_outflow = problem[:source_user] + F = problem[:F] + + for node_id in only(constraints_outflow.axes) + constraint = constraints_outflow[node_id] + user_idx = findsorted(user_demand.node_id, node_id) + capacity = + JuMP.normalized_rhs(constraint) + + user_demand.return_factor[user_idx] * + JuMP.value(F[(inflow_id(graph, node_id), node_id)]) + + JuMP.set_normalized_rhs(constraint, capacity) + end + + return nothing +end + """ Set the demand of the flow demand nodes. 2 cases: - Before the first allocation solve, set the demands to their full value; @@ -728,6 +770,8 @@ function allocate_priority!( priorities[priority_idx], collect_demands, ) + + adjust_capacities_returnflow!(allocation_model, p) return nothing end @@ -756,6 +800,8 @@ function allocate!( end end + set_initial_capacities_returnflow!(allocation_model) + for priority_idx in eachindex(priorities) allocate_priority!(allocation_model, u, p, t, priority_idx; collect_demands) end diff --git a/core/src/util.jl b/core/src/util.jl index 25dfbdee1..3b2131d1e 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -588,30 +588,6 @@ is_flow_constraining(node::AbstractParameterNode) = hasfield(typeof(node), :max_ is_flow_direction_constraining(node::AbstractParameterNode) = (nameof(typeof(node)) ∈ [:Pump, :Outlet, :TabulatedRatingCurve, :FractionalFlow]) -"""Find out whether a path exists between a start node and end node in the given allocation network.""" -function allocation_path_exists_in_graph( - graph::MetaGraph, - start_node_id::NodeID, - end_node_id::NodeID, -)::Bool - node_ids_visited = Set{NodeID}() - stack = [start_node_id] - - while !isempty(stack) - current_node_id = pop!(stack) - if current_node_id == end_node_id - return true - end - if !(current_node_id in node_ids_visited) - push!(node_ids_visited, current_node_id) - for outneighbor_node_id in outflow_ids_allocation(graph, current_node_id) - push!(stack, outneighbor_node_id) - end - end - end - return false -end - function has_main_network(allocation::Allocation)::Bool if !is_active(allocation) false diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 4d503e92e..996e4ed83 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -29,18 +29,19 @@ u = ComponentVector(; storage = zeros(length(p.basin.node_id))) Ribasim.allocate!(p, allocation_model, 0.0, u) + # Last priority (= 2) flows F = allocation_model.problem[:F] @test JuMP.value(F[(NodeID(:Basin, 2), NodeID(:Basin, 6))]) ≈ 0.0 @test JuMP.value(F[(NodeID(:Basin, 2), NodeID(:UserDemand, 10))]) ≈ 0.5 - @test JuMP.value(F[(NodeID(:Basin, 8), NodeID(:UserDemand, 12))]) ≈ 0.0 - @test JuMP.value(F[(NodeID(:Basin, 6), NodeID(:Basin, 8))]) ≈ 0.0 + @test JuMP.value(F[(NodeID(:Basin, 8), NodeID(:UserDemand, 12))]) ≈ 2.0 + @test JuMP.value(F[(NodeID(:Basin, 6), NodeID(:Basin, 8))]) ≈ 2.0 @test JuMP.value(F[(NodeID(:FlowBoundary, 1), NodeID(:Basin, 2))]) ≈ 0.5 @test JuMP.value(F[(NodeID(:Basin, 6), NodeID(:UserDemand, 11))]) ≈ 0.0 allocated = p.user_demand.allocated @test allocated[1] ≈ [0.0, 0.5] @test allocated[2] ≈ [4.0, 0.0] - @test allocated[3] ≈ [0.0, 0.0] + @test allocated[3] ≈ [0.0, 2.0] # Test getting and setting UserDemand demands (; user_demand) = p @@ -226,7 +227,7 @@ end @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))] ≈ - [0.001, 0.002, 0.002] + [0.001, 0.002, 0.0011] # Solving for the main network, containing subnetworks as UserDemands allocation_model = allocation_models[1] @@ -331,7 +332,7 @@ end @test ispath(toml_path) model = Ribasim.Model(toml_path) (; p) = model.integrator - (; graph, allocation, flow_demand) = p + (; graph, allocation, flow_demand, user_demand, level_boundary) = p # Test has_external_demand @test !any( @@ -373,14 +374,14 @@ end t = 0.0 - # Priority 1 + ## Priority 1 Ribasim.allocate_priority!(allocation_model, model.integrator.u, p, t, 1) objective = JuMP.objective_function(problem) @test F_abs_flow_demand[node_id_with_flow_demand] in keys(objective.terms) @test flow_demand.demand[1] == flow_demand.demand_itp[1](t) @test JuMP.normalized_rhs(constraint_flow_out) == Inf - # Priority 2 + ## Priority 2 Ribasim.allocate_priority!(allocation_model, model.integrator.u, p, t, 2) # The flow at priority 1 through the node with a flow demand at priority 2 # is subtracted from this flow demand @@ -389,7 +390,7 @@ end Ribasim.get_user_demand(p, NodeID(NodeType.UserDemand, 6), 1) @test JuMP.normalized_rhs(constraint_flow_out) == 0.0 - # Priority 3 + ## Priority 3 Ribasim.allocate_priority!(allocation_model, model.integrator.u, p, t, 3) @test JuMP.normalized_rhs(constraint_flow_out) == Inf # The flow from the source is used up in previous priorities @@ -397,6 +398,17 @@ end # So flow from the flow buffer is used for UserDemand #4 @test JuMP.value(F_flow_buffer_out[node_id_with_flow_demand]) == Ribasim.get_user_demand(p, NodeID(NodeType.UserDemand, 4), 3) + # Flow taken from buffer + @test JuMP.value(only(F_flow_buffer_out)) == user_demand.demand_itp[1][3](t) + # No flow coming from level boundary + @test JuMP.value(F[(only(level_boundary.node_id), node_id_with_flow_demand)]) == 0 + + ## Priority 4 + Ribasim.allocate_priority!(allocation_model, model.integrator.u, p, t, 4) + # Get demand from buffers + d = user_demand.demand_itp[3][4](t) + @assert JuMP.value(F[(NodeID(NodeType.UserDemand, 4), NodeID(NodeType.Basin, 7))]) + + JuMP.value(F[(NodeID(NodeType.UserDemand, 6), NodeID(NodeType.Basin, 7))]) == d end @testitem "flow_demand_with_max_flow_rate" begin diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index 92fa67ca1..dabf834f4 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -122,8 +122,9 @@ The capacities are determined in different ways: There are also capacities for special edges: -- $C^B_S \in \mathbb{R}^b_{\ge 0}$ where $b = \# B_S$ is the number of basins, for the flow supplied by basins. +- $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. # The optimization problem @@ -198,25 +199,23 @@ $$ When performing subnetwork demand collection, these capacities are set to $\infty$ for edges which connect the main network to a subnetwork. ::: -Similar constraints hold for the flow out of basins and flow demand buffers: +Similar constraints hold for the flow out of basins, flow demand buffers and user demand outflow sources: $$ -F^\text{basin out}_{i} \le (C^B_S)_i, \quad \forall i \in B_S, +F^\text{basin out}_{i} \le (C^{FD}_S)_i, \quad \forall i \in B_S, $$ $$ -F^\text{buffer out}_{i} \le (C^FD_S)_i, \quad \forall i \in FD_S. +F^\text{buffer out}_{i} \le (C^{FD}_S)_i, \quad \forall i \in FD_S, $$ - -- `UserDemand` outflow: The outflow of the UserDemand is dictated by the inflow and the return factor: $$ - F_{ik} = r_k \cdot F_{kj} \quad - \quad \forall k \in U_S, \quad - V^{\text{in}}_S(k) = \{i\},\; - V^{\text{out}}_S(k) = \{j\}. -$$ {#eq-returnflowconstraint} -Here we use that each UserDemand node in the allocation network has a unique in-edge and out-edge. -- User demand: UserDemand demand constraints are discussed in [the next section](allocation.qmd#sec-solving-allocation). +F_{ij} \le (C^{UD}_S)_i, \quad \forall i \in U_S, \quad V_S^{\text{out}}(i) = \{j\}. +$$ +Here we use that each UserDemand node in the allocation network has a unique outflow edge. The user outflow source capacities are increased after each optimization solve by the return fraction: +$$ + r_i \cdot F_{ki}, \quad V_S^{\text{in}}(i) = \{k\}. +$$ + - Fractional flow: Let $L_S \subset V_S$ be the set of nodes in the max flow graph with fractional flow outneighbors, and $f_j$ the flow fraction associated with fractional flow node $j \in V_S$. Then $$ F_{ij} \le f_j \sum_{k\in V^\text{in}_S(i)} F_{ki} \qquad @@ -226,37 +225,6 @@ $$ {#eq-fractionalflowconstraint} - Flow sign: Furthermore there are the non-negativity constraints for the flows and allocations, see [The optimization variables](allocation.qmd#the-optimization-variables). -## Final notes on the allocation problem - -### UserDemands using their own return flow - -If not explicitly avoided, UserDemands can use their own return flow in this allocation problem formulation. -Therefore, return flow of UserDemand is only taken into account by allocation if that return flow is downstream of the UserDemand where it comes from. That is, if there is no path in the directed allocation network from the UserDemand outflow node back to the UserDemand. - -# Solving the allocation problem {#sec-solving-allocation} - -The allocation problem for an allocation network at time $t$ is solved per priority, as follows: - -1. Define a capacity matrix with capacities as described above, that will be updated for each priority: -$$ - C_S^p \leftarrow C_S; -$$ -2. Set the capacities of the edges that end in a UserDemand to their priority 1 demands: -$$ - (C_S^p)_{i,j} \leftarrow d_j^1(t) \quad\text{ for all } (i,j) \in U_S; -$$ -3. Maximize the objective function given the constraints; -4. Subtract the used flows from the edge capacities: -$$ - C_S^p \leftarrow C_S^p - F; -$$ -5. Repeat steps 2-4 for the remaining priorities up to $p_{\max}$. - -:::{.callout-note} -In the future there will be 2 more optimization solves: -- One before optimizing for demands, taking the demand/supply from basins into account; -- One after optimizing for demands, taking preferences over sources into account. -::: ## Example diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index 243e0d9f2..e3da5d330 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -501,7 +501,7 @@ min_level | Float64 | $m$ | - max_level | Float64 | $m$ | - priority | Int32 | - | positive -# FlowDemand {#sec-level_demand} +# FlowDemand {#sec-flow_demand} A `FlowDemand` node associates a non-consuming flow demand to a connector node (e.g. `Pump`, `TabulatedRatingCurve`) for one single priority. FlowDemand nodes can set a flow demand only for a single connector node. diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 3509ead1f..2ff7bac8b 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -121,7 +121,7 @@ def subnetwork_model() -> Model: Node(11, Point(3, 3), subnetwork_id=2), [ user_demand.Static( - demand=[5.0], return_factor=0.9, min_level=0.9, priority=1 + demand=[5.0], return_factor=0.5, min_level=0.9, priority=1 ) ], ) @@ -939,7 +939,7 @@ def flow_demand_model() -> Model: Node(8, Point(3, -2), subnetwork_id=2), [ user_demand.Static( - priority=[4], demand=1e-3, return_factor=1.0, min_level=0.2 + priority=[4], demand=2e-3, return_factor=1.0, min_level=0.2 ) ], )