Skip to content

Commit

Permalink
Treat UserDemand return flow as a source (#1226)
Browse files Browse the repository at this point in the history 
Fixes #1212.
  • Loading branch information
@SouthEndMusic
SouthEndMusic authored Mar 21, 2024
1 parent 503443d commit 040afee
Show file tree
Hide file tree
Showing 7 changed files with 107 additions and 142 deletions.
89 changes: 26 additions & 63 deletions core/src/allocation_init.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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)
Expand Down
46 changes: 46 additions & 0 deletions core/src/allocation_optim.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -728,6 +770,8 @@ function allocate_priority!(
priorities[priority_idx],
collect_demands,
)

adjust_capacities_returnflow!(allocation_model, p)
return nothing
end

Expand Down Expand Up @@ -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
Expand Down
24 changes: 0 additions & 24 deletions core/src/util.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
28 changes: 20 additions & 8 deletions core/test/allocation_test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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
Expand All @@ -389,14 +390,25 @@ 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
@test JuMP.value(F[(NodeID(NodeType.LevelBoundary, 1), node_id_with_flow_demand)]) == 0
# 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
Expand Down
56 changes: 12 additions & 44 deletions docs/core/allocation.qmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
Loading

0 comments on commit 040afee

Please sign in to comment.