From 427eeb44b9332282c8c55266938be59847055cc1 Mon Sep 17 00:00:00 2001 From: Bart de Koning <74617371+SouthEndMusic@users.noreply.github.com> Date: Fri, 16 Feb 2024 12:57:14 +0100 Subject: [PATCH] Water supply/demand by basins in allocation (new node type `AllocationTarget`) (#1082) Fixes https://github.com/Deltares/Ribasim/issues/564. --- core/src/Ribasim.jl | 3 +- .../src/{allocation.jl => allocation_init.jl} | 621 +++++------------- core/src/allocation_optim.jl | 595 +++++++++++++++++ core/src/callback.jl | 6 +- core/src/parameter.jl | 25 +- core/src/read.jl | 37 +- core/src/schema.jl | 17 + core/src/solve.jl | 2 +- core/src/util.jl | 32 +- core/src/validation.jl | 6 +- core/test/allocation_test.jl | 100 ++- docs/core/allocation.qmd | 82 ++- docs/core/usage.qmd | 169 +++-- python/ribasim/ribasim/__init__.py | 2 + python/ribasim/ribasim/config.py | 13 + python/ribasim/ribasim/geometry/node.py | 2 + python/ribasim/ribasim/model.py | 2 + python/ribasim/ribasim/schemas.py | 15 + .../ribasim_testmodels/__init__.py | 2 + .../ribasim_testmodels/allocation.py | 123 +++- ribasim_qgis/core/nodes.py | 42 ++ 21 files changed, 1328 insertions(+), 568 deletions(-) rename core/src/{allocation.jl => allocation_init.jl} (59%) create mode 100644 core/src/allocation_optim.jl diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index 61684674a..5d4c6c8b6 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -82,7 +82,8 @@ include("parameter.jl") include("validation.jl") include("solve.jl") include("logging.jl") -include("allocation.jl") +include("allocation_init.jl") +include("allocation_optim.jl") include("util.jl") include("sparsity.jl") include("graph.jl") diff --git a/core/src/allocation.jl b/core/src/allocation_init.jl similarity index 59% rename from core/src/allocation.jl rename to core/src/allocation_init.jl index 568cd636f..39cfd100b 100644 --- a/core/src/allocation.jl +++ b/core/src/allocation_init.jl @@ -1,7 +1,7 @@ """Find the edges from the main network to a subnetwork.""" function find_subnetwork_connections!(p::Parameters)::Nothing - (; allocation, graph, user) = p - n_priorities = length(user.priorities) + (; allocation, graph, allocation) = p + n_priorities = length(allocation.priorities) (; subnetwork_demands, subnetwork_allocateds) = allocation for node_id in graph[].node_ids[1] for outflow_id in outflow_ids(graph, node_id) @@ -367,6 +367,27 @@ function add_variables_flow!( return nothing end +""" +Add the variables for supply/demand of a basin to the problem. +The variable indices are the node_ids of the basins in the subnetwork. +""" +function add_variables_basin!( + problem::JuMP.Model, + p::Parameters, + allocation_network_id::Int, +)::Nothing + (; graph) = p + node_ids_basin = [ + node_id for node_id in graph[].node_ids[allocation_network_id] if + graph[node_id].type == :basin + ] + problem[:F_basin_in] = + JuMP.@variable(problem, F_basin_in[node_id = node_ids_basin,] >= 0.0) + problem[:F_basin_out] = + JuMP.@variable(problem, F_basin_out[node_id = node_ids_basin,] >= 0.0) + return nothing +end + """ Certain allocation distribution types use absolute values in the objective function. Since most optimization packages do not support the absolute value function directly, @@ -383,7 +404,17 @@ function add_variables_absolute_value!( (; main_network_connections) = allocation if startswith(config.allocation.objective_type, "linear") node_ids = graph[].node_ids[allocation_network_id] - node_ids_user = [node_id for node_id in node_ids if node_id.type == NodeType.User] + node_ids_user = NodeID[] + node_ids_basin = NodeID[] + + for node_id in node_ids + type = node_id.type + if type == NodeType.User + push!(node_ids_user, node_id) + elseif type == NodeType.Basin + push!(node_ids_basin, node_id) + end + end # For the main network, connections to subnetworks are treated as users if is_main_network(allocation_network_id) @@ -394,7 +425,9 @@ function add_variables_absolute_value!( end end - problem[:F_abs] = JuMP.@variable(problem, F_abs[node_id = node_ids_user]) + problem[:F_abs_user] = JuMP.@variable(problem, F_abs_user[node_id = node_ids_user]) + problem[:F_abs_basin] = + JuMP.@variable(problem, F_abs_basin[node_id = node_ids_basin]) end return nothing end @@ -489,6 +522,30 @@ function outflow_ids_allocation(graph::MetaGraph, node_id::NodeID) return outflow_ids end +function get_basin_inflow( + problem::JuMP.Model, + node_id::NodeID, +)::Union{JuMP.VariableRef, Float64} + F_basin_in = problem[:F_basin_in] + return if node_id in only(F_basin_in.axes) + F_basin_in[node_id] + else + 0.0 + end +end + +function get_basin_outflow( + problem::JuMP.Model, + node_id::NodeID, +)::Union{JuMP.VariableRef, Float64} + F_basin_out = problem[:F_basin_out] + return if node_id in only(F_basin_out.axes) + F_basin_out[node_id] + else + 0.0 + end +end + """ Add the flow conservation constraints to the allocation problem. The constraint indices are user node IDs. @@ -514,10 +571,11 @@ function add_constraints_flow_conservation!( problem[:flow_conservation] = JuMP.@constraint( problem, [node_id = node_ids_conservation], - sum([ + get_basin_inflow(problem, node_id) + sum([ F[(node_id, outneighbor_id)] for outneighbor_id in outflow_ids_allocation(graph, node_id) - ]) == sum([ + ]) == + get_basin_outflow(problem, node_id) + sum([ F[(inneighbor_id, node_id)] for inneighbor_id in inflow_ids_allocation(graph, node_id) ]), @@ -564,71 +622,106 @@ expr_abs >= expr expr_abs >= -expr """ function add_constraints_absolute_value!( + problem::JuMP.Model, + flow_per_node::Dict{NodeID, JuMP.VariableRef}, + F_abs::JuMP.Containers.DenseAxisArray, + objective_type::String, + variable_type::String, +)::Nothing + # Example demand + d = 2.0 + + node_ids = only(F_abs.axes) + + if objective_type == "linear_absolute" + # These constraints together make sure that F_abs_* acts as the absolute + # value F_abs_* = |x| where x = F-d (here for example d = 2) + base_name = "abs_positive_$variable_type" + problem[Symbol(base_name)] = JuMP.@constraint( + problem, + [node_id = node_ids], + F_abs[node_id] >= (flow_per_node[node_id] - d), + base_name = base_name + ) + base_name = "abs_negative_$variable_type" + problem[Symbol(base_name)] = JuMP.@constraint( + problem, + [node_id = node_ids], + F_abs[node_id] >= -(flow_per_node[node_id] - d), + base_name = base_name + ) + elseif objective_type == "linear_relative" + # These constraints together make sure that F_abs_user acts as the absolute + # value F_abs_user = |x| where x = 1-F/d (here for example d = 2) + base_name = "abs_positive_$variable_type" + problem[Symbol(base_name)] = JuMP.@constraint( + problem, + [node_id = node_ids], + F_abs[node_id] >= (1 - flow_per_node[node_id] / d), + base_name = base_name + ) + base_name = "abs_negative_$variable_type" + problem[Symbol(base_name)] = JuMP.@constraint( + problem, + [node_id = node_ids], + F_abs[node_id] >= -(1 - flow_per_node[node_id] / d), + base_name = base_name + ) + end + return nothing +end + +""" +Add constraints so that variables F_abs_user act as the +absolute value of the expression comparing flow to an user to its demand. +""" +function add_constraints_absolute_value_user!( problem::JuMP.Model, p::Parameters, - allocation_network_id::Int, config::Config, )::Nothing - (; graph, allocation) = p - (; main_network_connections) = allocation + (; graph) = p objective_type = config.allocation.objective_type if startswith(objective_type, "linear") - node_ids = graph[].node_ids[allocation_network_id] - node_ids_user = [node_id for node_id in node_ids if node_id.type == NodeType.User] + F = problem[:F] + F_abs_user = problem[:F_abs_user] - # For the main network, connections to subnetworks are treated as users - if is_main_network(allocation_network_id) - for connections_subnetwork in main_network_connections - for connection in connections_subnetwork - push!(node_ids_user, connection[2]) - end - end - end + flow_per_node = Dict( + node_id => F[(only(inflow_ids_allocation(graph, node_id)), node_id)] for + node_id in only(F_abs_user.axes) + ) - node_ids_user_inflow = Dict( - node_id_user => only(inflow_ids_allocation(graph, node_id_user)) for - node_id_user in node_ids_user + add_constraints_absolute_value!( + problem, + flow_per_node, + F_abs_user, + objective_type, + "user", + ) + end + return nothing +end + +""" +Add constraints so that variables F_abs_basin act as the +absolute value of the expression comparing flow to a basin to its demand. +""" +function add_constraints_absolute_value_basin!(problem::JuMP.Model, config::Config)::Nothing + objective_type = config.allocation.objective_type + if startswith(objective_type, "linear") + F_basin_in = problem[:F_basin_in] + F_abs_basin = problem[:F_abs_basin] + flow_per_node = + Dict(node_id => F_basin_in[node_id] for node_id in only(F_abs_basin.axes)) + + add_constraints_absolute_value!( + problem, + flow_per_node, + F_abs_basin, + objective_type, + "basin", ) - F = problem[:F] - F_abs = problem[:F_abs] - d = 2.0 - - if config.allocation.objective_type == "linear_absolute" - # These constraints together make sure that F_abs acts as the absolute - # value F_abs = |x| where x = F-d (here for example d = 2) - problem[:abs_positive] = JuMP.@constraint( - problem, - [node_id_user = node_ids_user], - F_abs[node_id_user] >= - (F[(node_ids_user_inflow[node_id_user], node_id_user)] - d), - base_name = "abs_positive" - ) - problem[:abs_negative] = JuMP.@constraint( - problem, - [node_id_user = node_ids_user], - F_abs[node_id_user] >= - -(F[(node_ids_user_inflow[node_id_user], node_id_user)] - d), - base_name = "abs_negative" - ) - elseif config.allocation.objective_type == "linear_relative" - # These constraints together make sure that F_abs acts as the absolute - # value F_abs = |x| where x = 1-F/d (here for example d = 2) - problem[:abs_positive] = JuMP.@constraint( - problem, - [node_id_user = node_ids_user], - F_abs[node_id_user] >= - (1 - F[(node_ids_user_inflow[node_id_user], node_id_user)] / d), - base_name = "abs_positive" - ) - problem[:abs_negative] = JuMP.@constraint( - problem, - [node_id_user = node_ids_user], - F_abs[node_id_user] >= - -(1 - F[(node_ids_user_inflow[node_id_user], node_id_user)] / d), - base_name = "abs_negative" - ) - end end return nothing end @@ -684,6 +777,24 @@ function add_constraints_fractional_flow!( return nothing end +""" +Add the basin flow constraints to the allocation problem. +The constraint indices are the basin node IDs. + +Constraint: +flow out of basin <= basin capacity +""" +function add_constraints_basin_flow!(problem::JuMP.Model)::Nothing + F_basin_out = problem[:F_basin_out] + problem[:basin_outflow] = JuMP.@constraint( + problem, + [node_id = only(F_basin_out.axes)], + F_basin_out[node_id] <= 0.0, + base_name = "basin_outflow" + ) + return nothing +end + """ Construct the allocation problem for the current subnetwork as a JuMP.jl model. """ @@ -698,16 +809,18 @@ function allocation_problem( # Add variables to problem add_variables_flow!(problem, p, allocation_network_id) + add_variables_basin!(problem, p, allocation_network_id) add_variables_absolute_value!(problem, p, allocation_network_id, config) - # TODO: Add variables for allocation to basins # Add constraints to problem add_constraints_capacity!(problem, capacity, p, allocation_network_id) add_constraints_source!(problem, p, allocation_network_id) add_constraints_flow_conservation!(problem, p, allocation_network_id) add_constraints_user_returnflow!(problem, p, allocation_network_id) - add_constraints_absolute_value!(problem, p, allocation_network_id, config) + add_constraints_absolute_value_user!(problem, p, config) + add_constraints_absolute_value_basin!(problem, config) add_constraints_fractional_flow!(problem, p, allocation_network_id) + add_constraints_basin_flow!(problem) return problem end @@ -745,385 +858,3 @@ function AllocationModel( Δt_allocation, ) end - -""" -Add a term to the expression of the objective function corresponding to -the demand of a user. -""" -function add_user_term!( - ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, - edge::Tuple{NodeID, NodeID}, - objective_type::Symbol, - demand::Float64, - model::AllocationModel, -)::Nothing - (; problem) = model - F = problem[:F] - F_edge = F[edge] - node_id_user = edge[2] - - if objective_type == :quadratic_absolute - # Objective function ∑ (F - d)^2 - JuMP.add_to_expression!(ex, 1, F_edge, F_edge) - JuMP.add_to_expression!(ex, -2 * demand, F_edge) - JuMP.add_to_expression!(ex, demand^2) - - elseif objective_type == :quadratic_relative - # Objective function ∑ (1 - F/d)^2 - if demand ≈ 0 - return nothing - end - JuMP.add_to_expression!(ex, 1.0 / demand^2, F_edge, F_edge) - JuMP.add_to_expression!(ex, -2.0 / demand, F_edge) - JuMP.add_to_expression!(ex, 1.0) - - elseif objective_type == :linear_absolute - # Objective function ∑ |F - d| - JuMP.set_normalized_rhs(problem[:abs_positive][node_id_user], -demand) - JuMP.set_normalized_rhs(problem[:abs_negative][node_id_user], demand) - - elseif objective_type == :linear_relative - # Objective function ∑ |1 - F/d| - JuMP.set_normalized_coefficient( - problem[:abs_positive][node_id_user], - F_edge, - iszero(demand) ? 0 : 1 / demand, - ) - JuMP.set_normalized_coefficient( - problem[:abs_negative][node_id_user], - F_edge, - iszero(demand) ? 0 : -1 / demand, - ) - else - error("Invalid allocation objective type $objective_type.") - end - return nothing -end - -""" -Set the objective for the given priority. -For an objective with absolute values this also involves adjusting constraints. -""" -function set_objective_priority!( - allocation_model::AllocationModel, - p::Parameters, - t::Float64, - priority_idx::Int, -)::Nothing - (; objective_type, problem, allocation_network_id) = allocation_model - (; graph, user, allocation) = p - (; demand, demand_itp, demand_from_timeseries, node_id) = user - (; main_network_connections, subnetwork_demands) = allocation - edge_ids = graph[].edge_ids[allocation_network_id] - - F = problem[:F] - if objective_type in [:quadratic_absolute, :quadratic_relative] - ex = JuMP.QuadExpr() - elseif objective_type in [:linear_absolute, :linear_relative] - ex = sum(problem[:F_abs]) - end - - # Terms for subnetworks as users - if is_main_network(allocation_network_id) - for connections_subnetwork in main_network_connections - for connection in connections_subnetwork - d = subnetwork_demands[connection][priority_idx] - add_user_term!(ex, connection, objective_type, d, allocation_model) - end - end - end - - # Terms for user nodes - for edge_id in edge_ids - node_id_user = edge_id[2] - if node_id_user.type != NodeType.User - continue - end - - user_idx = findsorted(node_id, node_id_user) - - if demand_from_timeseries[user_idx] - d = demand_itp[user_idx][priority_idx](t) - set_user_demand!(user, node_id_user, priority_idx, d) - else - d = get_user_demand(user, node_id_user, priority_idx) - end - - add_user_term!(ex, edge_id, objective_type, d, allocation_model) - end - - new_objective = JuMP.@expression(problem, ex) - JuMP.@objective(problem, Min, new_objective) - return nothing -end - -""" -Assign the allocations to the users as determined by the solution of the allocation problem. -""" -function assign_allocations!( - allocation_model::AllocationModel, - p::Parameters, - t::Float64, - priority_idx::Int; - collect_demands::Bool = false, -)::Nothing - (; problem, allocation_network_id) = allocation_model - (; graph, user, allocation) = p - (; - subnetwork_demands, - subnetwork_allocateds, - allocation_network_ids, - main_network_connections, - ) = allocation - (; record) = user - edge_ids = graph[].edge_ids[allocation_network_id] - main_network_source_edges = get_main_network_connections(p, allocation_network_id) - F = problem[:F] - for edge_id in edge_ids - # If this edge is a source edge from the main network to a subnetwork, - # and demands are being collected, add its flow to the demand of this edge - if collect_demands && - graph[edge_id...].allocation_network_id_source == allocation_network_id && - edge_id ∈ main_network_source_edges - allocated = JuMP.value(F[edge_id]) - subnetwork_demands[edge_id][priority_idx] += allocated - end - - user_node_id = edge_id[2] - - if user_node_id.type == NodeType.User - allocated = JuMP.value(F[edge_id]) - user_idx = findsorted(user.node_id, user_node_id) - user.allocated[user_idx][priority_idx] = allocated - - # Save allocations to record - push!(record.time, t) - push!(record.subnetwork_id, allocation_model.allocation_network_id) - push!(record.user_node_id, Int(user_node_id)) - push!(record.priority, user.priorities[priority_idx]) - push!(record.demand, user.demand[user_idx]) - push!(record.allocated, allocated) - # TODO: This is now the last abstraction before the allocation update, - # should be the average abstraction since the last allocation solve - push!( - record.abstracted, - get_flow(graph, inflow_id(graph, user_node_id), user_node_id, 0), - ) - end - end - - # Write the flows to the subnetworks as allocated flows - # in the allocation object - if is_main_network(allocation_network_id) - for (allocation_network_id, main_network_source_edges) in - zip(allocation_network_ids, main_network_connections) - if is_main_network(allocation_network_id) - continue - end - for edge_id in main_network_source_edges - subnetwork_allocateds[edge_id][priority_idx] = JuMP.value(F[edge_id]) - end - end - end - return nothing -end - -""" -Adjust the source flows. -""" -function adjust_source_capacities!( - allocation_model::AllocationModel, - p::Parameters, - priority_idx::Int; - collect_demands::Bool = false, -)::Nothing - (; problem) = allocation_model - (; graph, allocation) = p - (; allocation_network_id) = allocation_model - (; subnetwork_allocateds) = allocation - edge_ids = graph[].edge_ids[allocation_network_id] - source_constraints = problem[:source] - F = problem[:F] - - main_network_source_edges = get_main_network_connections(p, allocation_network_id) - - for edge_id in edge_ids - if graph[edge_id...].allocation_network_id_source == allocation_network_id - # If it is a source edge for this allocation problem - if priority_idx == 1 - # If the optimization was just started, i.e. sources have to be reset - if edge_id in main_network_source_edges - if collect_demands - # Set the source capacity to effectively unlimited if subnetwork demands are being collected - source_capacity = Inf - else - # Set the source capacity to the value allocated to the subnetwork over this edge - source_capacity = subnetwork_allocateds[edge_id][priority_idx] - end - else - # Reset the source to the current flow from the physical layer. - source_capacity = get_flow(graph, edge_id..., 0) - end - JuMP.set_normalized_rhs( - source_constraints[edge_id], - # It is assumed that the allocation procedure does not have to be differentiated. - source_capacity, - ) - else - # Subtract the allocated flow from the source. - JuMP.set_normalized_rhs( - source_constraints[edge_id], - JuMP.normalized_rhs(source_constraints[edge_id]) - - JuMP.value(F[edge_id]), - ) - end - end - end - return nothing -end - -""" -Set the values of the edge capacities. 2 cases: -- Before the first allocation solve, set the edge capacities to their full capacity; -- Before an allocation solve, subtract the flow used by allocation for the previous priority - from the edge capacities. -""" -function adjust_edge_capacities!( - allocation_model::AllocationModel, - p::Parameters, - priority_idx::Int, -)::Nothing - (; graph) = p - (; problem, capacity, allocation_network_id) = allocation_model - edge_ids = graph[].edge_ids[allocation_network_id] - constraints_capacity = problem[:capacity] - F = problem[:F] - - main_network_source_edges = get_main_network_connections(p, allocation_network_id) - - for edge_id in edge_ids - c = capacity[edge_id...] - - # These edges have no capacity constraints: - # - With infinite capacity - # - Being a source from the main network to a subnetwork - if isinf(c) || edge_id ∈ main_network_source_edges - continue - end - - if priority_idx == 1 - # Before the first allocation solve, set the edge capacities to their full capacity - JuMP.set_normalized_rhs(constraints_capacity[edge_id], c) - else - # Before an allocation solve, subtract the flow used by allocation for the previous priority - # from the edge capacities - JuMP.set_normalized_rhs( - constraints_capacity[edge_id], - JuMP.normalized_rhs(constraints_capacity[edge_id]) - JuMP.value(F[edge_id]), - ) - end - end -end - -""" -Save the allocation flows per physical edge. -""" -function save_allocation_flows!( - p::Parameters, - t::Float64, - allocation_model::AllocationModel, - priority::Int, - collect_demands::Bool, -)::Nothing - (; problem, allocation_network_id) = allocation_model - (; allocation, graph) = p - (; record) = allocation - F = problem[:F] - - for allocation_edge in first(F.axes) - flow = JuMP.value(F[allocation_edge]) - edge_metadata = graph[allocation_edge...] - (; node_ids) = edge_metadata - - for i in eachindex(node_ids)[1:(end - 1)] - push!(record.time, t) - push!(record.edge_id, edge_metadata.id) - push!(record.from_node_id, node_ids[i]) - push!(record.to_node_id, node_ids[i + 1]) - push!(record.subnetwork_id, allocation_network_id) - push!(record.priority, priority) - push!(record.flow, flow) - push!(record.collect_demands, collect_demands) - 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 users. -""" -function allocate!( - p::Parameters, - allocation_model::AllocationModel, - t::Float64; - collect_demands::Bool = false, -)::Nothing - (; user, allocation) = p - (; problem, allocation_network_id) = allocation_model - (; priorities) = user - (; subnetwork_demands) = allocation - - # TODO: Compute basin flow from vertical fluxes and basin volume. - # Set as basin demand if the net flow is negative, set as source - # in the flow_conservation constraints if the net flow is positive. - # Solve this as a separate problem before the priorities below - - main_network_source_edges = get_main_network_connections(p, allocation_network_id) - - if collect_demands - for main_network_connection in keys(subnetwork_demands) - if main_network_connection in main_network_source_edges - subnetwork_demands[main_network_connection] .= 0.0 - end - end - end - - for priority_idx in eachindex(priorities) - adjust_source_capacities!(allocation_model, p, priority_idx; collect_demands) - - # Subtract the flows used by the allocation of the previous priority from the capacities of the edges - # or set edge capacities if priority_idx = 1 - adjust_edge_capacities!(allocation_model, p, priority_idx) - - # Set the objective depending on the demands - # A new objective function is set instead of modifying the coefficients - # of an existing objective function because this is not supported for - # quadratic terms: - # https://jump.dev/JuMP.jl/v1.16/manual/objective/#Modify-an-objective-coefficient - set_objective_priority!(allocation_model, p, t, priority_idx) - - # Solve the allocation problem for this priority - JuMP.optimize!(problem) - @debug JuMP.solution_summary(problem) - if JuMP.termination_status(problem) !== JuMP.OPTIMAL - (; allocation_network_id) = allocation_model - priority = priorities[priority_idx] - error( - "Allocation of subnetwork $allocation_network_id, priority $priority coudn't find optimal solution.", - ) - end - - # Assign the allocations to the users for this priority - assign_allocations!(allocation_model, p, t, priority_idx; collect_demands) - - # Save the flows over all edges in the subnetwork - save_allocation_flows!( - p, - t, - allocation_model, - priorities[priority_idx], - collect_demands, - ) - end -end diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl new file mode 100644 index 000000000..812875639 --- /dev/null +++ b/core/src/allocation_optim.jl @@ -0,0 +1,595 @@ +""" +Add a term to the objective function given by the objective type, +depending in the provided flow variable and the associated demand. +""" +function add_objective_term!( + ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, + F_variable::JuMP.VariableRef, + demand::Float64, + objective_type::Symbol; + constraint_abs_positive::Union{JuMP.ConstraintRef, Nothing} = nothing, + constraint_abs_negative::Union{JuMP.ConstraintRef, Nothing} = nothing, +)::Nothing + if objective_type == :quadratic_absolute + # Objective function ∑ (F - d)^2 + JuMP.add_to_expression!(ex, 1, F_variable, F_variable) + JuMP.add_to_expression!(ex, -2 * demand, F_variable) + JuMP.add_to_expression!(ex, demand^2) + + elseif objective_type == :quadratic_relative + # Objective function ∑ (1 - F/d)^2 + if demand ≈ 0 + return nothing + end + JuMP.add_to_expression!(ex, 1.0 / demand^2, F_variable, F_variable) + JuMP.add_to_expression!(ex, -2.0 / demand, F_variable) + JuMP.add_to_expression!(ex, 1.0) + + elseif objective_type == :linear_absolute + # Objective function ∑ |F - d| + JuMP.set_normalized_rhs(constraint_abs_positive, -demand) + JuMP.set_normalized_rhs(constraint_abs_negative, demand) + + elseif objective_type == :linear_relative + # Objective function ∑ |1 - F/d| + JuMP.set_normalized_coefficient( + constraint_abs_positive, + F_variable, + iszero(demand) ? 0 : 1 / demand, + ) + JuMP.set_normalized_coefficient( + constraint_abs_negative, + F_variable, + iszero(demand) ? 0 : -1 / demand, + ) + else + error("Invalid allocation objective type $objective_type.") + end + return nothing +end + +""" +Add a term to the expression of the objective function corresponding to +the demand of a user. +""" +function add_user_term!( + ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, + edge::Tuple{NodeID, NodeID}, + objective_type::Symbol, + demand::Float64, + problem::JuMP.Model, +)::Nothing + F = problem[:F] + F_edge = F[edge] + node_id_user = edge[2] + + if objective_type in [:linear_absolute, :linear_relative] + constraint_abs_positive = problem[:abs_positive_user][node_id_user] + constraint_abs_negative = problem[:abs_negative_user][node_id_user] + else + constraint_abs_positive = nothing + constraint_abs_negative = nothing + end + + add_objective_term!( + ex, + F_edge, + demand, + objective_type; + constraint_abs_positive, + constraint_abs_negative, + ) +end + +""" +Add a term to the expression of the objective function corresponding to +the demand of a basin. +""" +function add_basin_term!( + ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, + problem::JuMP.Model, + demand::Float64, + objective_type::Symbol, + node_id::NodeID, +)::Nothing + F_basin_in = problem[:F_basin_in] + F_basin = F_basin_in[node_id] + + if objective_type in [:linear_absolute, :linear_relative] + constraint_abs_positive = problem[:abs_positive_basin][node_id] + constraint_abs_negative = problem[:abs_negative_basin][node_id] + else + constraint_abs_positive = nothing + constraint_abs_negative = nothing + end + + add_objective_term!( + ex, + F_basin, + demand, + objective_type; + constraint_abs_positive, + constraint_abs_negative, + ) + return nothing +end + +""" +Set the objective for the given priority. +For an objective with absolute values this also involves adjusting constraints. +""" +function set_objective_priority!( + allocation_model::AllocationModel, + p::Parameters, + u::ComponentVector, + t::Float64, + priority_idx::Int, +)::Nothing + (; objective_type, problem, allocation_network_id) = allocation_model + (; graph, user, allocation) = p + (; demand_itp, demand_from_timeseries, node_id) = user + (; main_network_connections, subnetwork_demands, priorities) = allocation + edge_ids = graph[].edge_ids[allocation_network_id] + + F = problem[:F] + if objective_type in [:quadratic_absolute, :quadratic_relative] + ex = JuMP.QuadExpr() + elseif objective_type in [:linear_absolute, :linear_relative] + ex = JuMP.AffExpr() + ex += sum(problem[:F_abs_user]) + ex += sum(problem[:F_abs_basin]) + end + + demand_max = 0.0 + + # Terms for subnetworks as users + if is_main_network(allocation_network_id) + for connections_subnetwork in main_network_connections + for connection in connections_subnetwork + d = subnetwork_demands[connection][priority_idx] + demand_max = max(demand_max, d) + add_user_term!(ex, connection, objective_type, d, problem) + end + end + end + + # Terms for user nodes + for edge_id in edge_ids + node_id_user = edge_id[2] + if node_id_user.type != NodeType.User + continue + end + + user_idx = findsorted(node_id, node_id_user) + if demand_from_timeseries[user_idx] + d = demand_itp[user_idx][priority_idx](t) + set_user_demand!(p, node_id_user, priority_idx, d) + else + d = get_user_demand(p, node_id_user, priority_idx) + end + demand_max = max(demand_max, d) + add_user_term!(ex, edge_id, objective_type, d, problem) + end + + # Terms for basins + F_basin_in = problem[:F_basin_in] + for node_id in only(F_basin_in.axes) + priority_basin = get_basin_priority(p, node_id) + priority_now = priorities[priority_idx] + d = + priority_basin == priority_now ? + get_basin_demand(allocation_model, u, p, t, node_id) : 0.0 + add_basin_term!(ex, problem, d, objective_type, node_id) + end + + new_objective = JuMP.@expression(problem, ex) + JuMP.@objective(problem, Min, new_objective) + return nothing +end + +""" +Assign the allocations to the users as determined by the solution of the allocation problem. +""" +function assign_allocations!( + allocation_model::AllocationModel, + p::Parameters, + t::Float64, + priority_idx::Int; + collect_demands::Bool = false, +)::Nothing + (; problem, allocation_network_id) = allocation_model + (; graph, user, allocation) = p + (; + subnetwork_demands, + subnetwork_allocateds, + allocation_network_ids, + main_network_connections, + ) = allocation + (; record) = user + edge_ids = graph[].edge_ids[allocation_network_id] + main_network_source_edges = get_main_network_connections(p, allocation_network_id) + F = problem[:F] + for edge_id in edge_ids + # If this edge is a source edge from the main network to a subnetwork, + # and demands are being collected, add its flow to the demand of this edge + if collect_demands && + graph[edge_id...].allocation_network_id_source == allocation_network_id && + edge_id ∈ main_network_source_edges + allocated = JuMP.value(F[edge_id]) + subnetwork_demands[edge_id][priority_idx] += allocated + end + + user_node_id = edge_id[2] + + if user_node_id.type == NodeType.User + allocated = JuMP.value(F[edge_id]) + user_idx = findsorted(user.node_id, user_node_id) + user.allocated[user_idx][priority_idx] = allocated + + # Save allocations to record + push!(record.time, t) + push!(record.subnetwork_id, allocation_model.allocation_network_id) + push!(record.user_node_id, Int(user_node_id)) + push!(record.priority, allocation.priorities[priority_idx]) + push!(record.demand, user.demand[user_idx]) + push!(record.allocated, allocated) + # TODO: This is now the last abstraction before the allocation update, + # should be the average abstraction since the last allocation solve + push!( + record.abstracted, + get_flow(graph, inflow_id(graph, user_node_id), user_node_id, 0), + ) + end + end + + # Write the flows to the subnetworks as allocated flows + # in the allocation object + if is_main_network(allocation_network_id) + for (allocation_network_id, main_network_source_edges) in + zip(allocation_network_ids, main_network_connections) + if is_main_network(allocation_network_id) + continue + end + for edge_id in main_network_source_edges + subnetwork_allocateds[edge_id][priority_idx] = JuMP.value(F[edge_id]) + end + end + end + return nothing +end + +""" +Adjust the source flows. +""" +function adjust_source_capacities!( + allocation_model::AllocationModel, + p::Parameters, + priority_idx::Int; + collect_demands::Bool = false, +)::Nothing + (; problem) = allocation_model + (; graph, allocation) = p + (; allocation_network_id) = allocation_model + (; subnetwork_allocateds) = allocation + edge_ids = graph[].edge_ids[allocation_network_id] + source_constraints = problem[:source] + F = problem[:F] + + main_network_source_edges = get_main_network_connections(p, allocation_network_id) + + for edge_id in edge_ids + if graph[edge_id...].allocation_network_id_source == allocation_network_id + # If it is a source edge for this allocation problem + if priority_idx == 1 + # If the optimization was just started, i.e. sources have to be reset + if edge_id in main_network_source_edges + if collect_demands + # Set the source capacity to effectively unlimited if subnetwork demands are being collected + source_capacity = Inf + else + # Set the source capacity to the value allocated to the subnetwork over this edge + source_capacity = subnetwork_allocateds[edge_id][priority_idx] + end + else + # Reset the source to the current flow from the physical layer. + source_capacity = get_flow(graph, edge_id..., 0) + end + JuMP.set_normalized_rhs( + source_constraints[edge_id], + # It is assumed that the allocation procedure does not have to be differentiated. + source_capacity, + ) + else + # Subtract the allocated flow from the source. + JuMP.set_normalized_rhs( + source_constraints[edge_id], + JuMP.normalized_rhs(source_constraints[edge_id]) - + JuMP.value(F[edge_id]), + ) + end + end + end + return nothing +end + +""" +Set the values of the edge capacities. 2 cases: +- Before the first allocation solve, set the edge capacities to their full capacity; +- Before an allocation solve, subtract the flow used by allocation for the previous priority + from the edge capacities. +""" +function adjust_edge_capacities!( + allocation_model::AllocationModel, + p::Parameters, + priority_idx::Int, +)::Nothing + (; graph) = p + (; problem, capacity, allocation_network_id) = allocation_model + edge_ids = graph[].edge_ids[allocation_network_id] + constraints_capacity = problem[:capacity] + F = problem[:F] + + main_network_source_edges = get_main_network_connections(p, allocation_network_id) + + for edge_id in edge_ids + c = capacity[edge_id...] + + # These edges have no capacity constraints: + # - With infinite capacity + # - Being a source from the main network to a subnetwork + if isinf(c) || edge_id ∈ main_network_source_edges + continue + end + + if priority_idx == 1 + # Before the first allocation solve, set the edge capacities to their full capacity + JuMP.set_normalized_rhs(constraints_capacity[edge_id], c) + else + # Before an allocation solve, subtract the flow used by allocation for the previous priority + # from the edge capacities + JuMP.set_normalized_rhs( + constraints_capacity[edge_id], + JuMP.normalized_rhs(constraints_capacity[edge_id]) - JuMP.value(F[edge_id]), + ) + end + end +end + +""" +Get several variables associated with a basin: +- Its current storage +- The allocation update interval +- The influx (sum of instantaneous vertical fluxes of the basin) +- The index of the connected allocation_target node (0 if such a + node does not exist) +- The index of the basin +""" +function get_basin_data( + allocation_model::AllocationModel, + p::Parameters, + u::ComponentVector, + node_id::NodeID, +) + (; graph, basin, allocation_target) = p + (; Δt_allocation) = allocation_model + @assert node_id.type == NodeType.Basin + influx = get_flow(graph, node_id, 0.0) + _, basin_idx = id_index(basin.node_id, node_id) + storage_basin = u.storage[basin_idx] + control_inneighbors = inneighbor_labels_type(graph, node_id, EdgeType.control) + if isempty(control_inneighbors) + allocation_target_idx = 0 + else + allocation_target_node_id = first(control_inneighbors) + allocation_target_idx = + findsorted(allocation_target.node_id, allocation_target_node_id) + end + return storage_basin, Δt_allocation, influx, allocation_target_idx, basin_idx +end + +""" +Get the capacity of the basin, i.e. the maximum +flow that can be abstracted from the basin if it is in a +state of surplus storage (0 if no reference levels are provided by +a allocation_target node). +Storages are converted to flows by dividing by the allocation timestep. +""" +function get_basin_capacity( + allocation_model::AllocationModel, + u::ComponentVector, + p::Parameters, + t::Float64, + node_id::NodeID, +)::Float64 + (; allocation_target) = p + @assert node_id.type == NodeType.Basin + storage_basin, Δt_allocation, influx, allocation_target_idx, basin_idx = + get_basin_data(allocation_model, p, u, node_id) + if iszero(allocation_target_idx) + return 0.0 + else + level_max = allocation_target.max_level[allocation_target_idx](t) + storage_max = get_storage_from_level(p.basin, basin_idx, level_max) + return max(0.0, (storage_basin - storage_max) / Δt_allocation + influx) + end +end + +""" +Get the demand of the basin, i.e. how large a flow the +basin needs to get to its minimum target level (0 if no +reference levels are provided by a allocation_target node). +Storages are converted to flows by dividing by the allocation timestep. +""" +function get_basin_demand( + allocation_model::AllocationModel, + u::ComponentVector, + p::Parameters, + t::Float64, + node_id::NodeID, +)::Float64 + (; allocation_target) = p + @assert node_id.type == NodeType.Basin + storage_basin, Δt_allocation, influx, allocation_target_idx, basin_idx = + get_basin_data(allocation_model, p, u, node_id) + if iszero(allocation_target_idx) + return 0.0 + else + level_min = allocation_target.min_level[allocation_target_idx](t) + storage_min = get_storage_from_level(p.basin, basin_idx, level_min) + return max(0.0, (storage_min - storage_basin) / Δt_allocation - influx) + end +end + +""" +Set the values of the basin outflows. 2 cases: +- Before the first allocation solve, set the capacities to their full capacity if there is surplus storage; +- Before an allocation solve, subtract the flow used by allocation for the previous priority + from the capacities. +""" +function adjust_basin_capacities!( + allocation_model::AllocationModel, + u::ComponentVector, + p::Parameters, + t::Float64, + priority_idx::Int, +)::Nothing + (; problem) = allocation_model + constraints_outflow = problem[:basin_outflow] + F_basin_out = problem[:F_basin_out] + + for node_id in only(constraints_outflow.axes) + constraint = constraints_outflow[node_id] + if priority_idx == 1 + JuMP.set_normalized_rhs( + constraint, + get_basin_capacity(allocation_model, u, p, t, node_id), + ) + else + JuMP.set_normalized_rhs( + constraint, + JuMP.normalized_rhs(constraint) - JuMP.value(F_basin_out[node_id]), + ) + end + end + + return nothing +end + +""" +Save the allocation flows per basin physical edge. +""" +function save_allocation_flows!( + p::Parameters, + t::Float64, + allocation_model::AllocationModel, + priority::Int, + collect_demands::Bool, +)::Nothing + (; problem, allocation_network_id) = allocation_model + (; allocation, graph) = p + (; record) = allocation + F = problem[:F] + F_basin_in = problem[:F_basin_in] + F_basin_out = problem[:F_basin_out] + + # Edge flows + for allocation_edge in first(F.axes) + flow = JuMP.value(F[allocation_edge]) + edge_metadata = graph[allocation_edge...] + (; node_ids) = edge_metadata + + for i in eachindex(node_ids)[1:(end - 1)] + push!(record.time, t) + push!(record.edge_id, edge_metadata.id) + push!(record.from_node_id, node_ids[i]) + push!(record.to_node_id, node_ids[i + 1]) + push!(record.subnetwork_id, allocation_network_id) + push!(record.priority, priority) + push!(record.flow, flow) + push!(record.collect_demands, collect_demands) + end + end + + # Basin flows + for node_id in graph[].node_ids[allocation_network_id] + if node_id.type == NodeType.Basin + flow = JuMP.value(F_basin_out[node_id]) - JuMP.value(F_basin_in[node_id]) + push!(record.time, t) + push!(record.edge_id, 0) + push!(record.from_node_id, node_id) + push!(record.to_node_id, node_id) + push!(record.subnetwork_id, allocation_network_id) + push!(record.priority, priority) + push!(record.flow, flow) + push!(record.collect_demands, collect_demands) + 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 users. +""" +function allocate!( + p::Parameters, + allocation_model::AllocationModel, + t::Float64, + u::ComponentVector; + collect_demands::Bool = false, +)::Nothing + (; allocation) = p + (; problem, allocation_network_id) = allocation_model + (; priorities, subnetwork_demands) = allocation + + main_network_source_edges = get_main_network_connections(p, allocation_network_id) + + if collect_demands + for main_network_connection in keys(subnetwork_demands) + if main_network_connection in main_network_source_edges + subnetwork_demands[main_network_connection] .= 0.0 + end + end + end + + for priority_idx in eachindex(priorities) + adjust_source_capacities!(allocation_model, p, priority_idx; collect_demands) + + # Subtract the flows used by the allocation of the previous priority from the capacities of the edges + # or set edge capacities if priority_idx = 1 + adjust_edge_capacities!(allocation_model, p, priority_idx) + + adjust_basin_capacities!(allocation_model, u, p, t, priority_idx) + + # Set the objective depending on the demands + # A new objective function is set instead of modifying the coefficients + # of an existing objective function because this is not supported for + # quadratic terms: + # https://jump.dev/JuMP.jl/v1.16/manual/objective/#Modify-an-objective-coefficient + set_objective_priority!(allocation_model, p, u, t, priority_idx) + + # Solve the allocation problem for this priority + JuMP.optimize!(problem) + @debug JuMP.solution_summary(problem) + if JuMP.termination_status(problem) !== JuMP.OPTIMAL + (; allocation_network_id) = allocation_model + priority = priorities[priority_idx] + error( + "Allocation of subnetwork $allocation_network_id, priority $priority coudn't find optimal solution.", + ) + end + + # Assign the allocations to the users for this priority + assign_allocations!(allocation_model, p, t, priority_idx; collect_demands) + + # Save the flows over all edges in the subnetwork + save_allocation_flows!( + p, + t, + allocation_model, + priorities[priority_idx], + collect_demands, + ) + end +end diff --git a/core/src/callback.jl b/core/src/callback.jl index 14c7ffafd..087f3e049 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -426,14 +426,14 @@ end "Solve the allocation problem for all users and assign allocated abstractions to user nodes." function update_allocation!(integrator)::Nothing - (; p, t) = integrator + (; p, t, u) = integrator (; allocation) = p (; allocation_models) = allocation # If a main network is present, collect demands of subnetworks if has_main_network(allocation) for allocation_model in Iterators.drop(allocation_models, 1) - allocate!(p, allocation_model, t; collect_demands = true) + allocate!(p, allocation_model, t, u; collect_demands = true) end end @@ -441,7 +441,7 @@ function update_allocation!(integrator)::Nothing # If a main network is present this is solved first, # which provides allocation to the subnetworks for allocation_model in allocation_models - allocate!(p, allocation_model, t) + allocate!(p, allocation_model, t, u) end end diff --git a/core/src/parameter.jl b/core/src/parameter.jl index d2306af19..7b0de5f93 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -63,6 +63,7 @@ allocation_network_ids: The unique sorted allocation network IDs allocation models: The allocation models for the main network and subnetworks corresponding to allocation_network_ids main_network_connections: (from_id, to_id) from the main network to the subnetwork per subnetwork +priorities: All used priority values. subnetwork_demands: The demand of an edge from the main network to a subnetwork record: A record of all flows computed by allocation optimization, eventually saved to output file @@ -71,6 +72,7 @@ struct Allocation allocation_network_ids::Vector{Int} allocation_models::Vector{AllocationModel} main_network_connections::Vector{Vector{Tuple{NodeID, NodeID}}} + priorities::Vector{Int} subnetwork_demands::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} subnetwork_allocateds::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} record::@NamedTuple{ @@ -456,13 +458,13 @@ struct PidControl{T} <: AbstractParameterNode end """ -demand: water flux demand of user per priority over time +demand: water flux demand of user per priority over time. + Each user has a demand for all priorities, + which is 0.0 if it is not provided explicitly. active: whether this node is active and thus demands water allocated: water flux currently allocated to user per priority return_factor: the factor in [0,1] of how much of the abstracted water is given back to the system min_level: The level of the source basin below which the user does not abstract -priorities: All used priority values. Each user has a demand for all these priorities, - which is 0.0 if it is not provided explicitly. record: Collected data of allocation optimizations for output file. """ struct User <: AbstractParameterNode @@ -474,7 +476,6 @@ struct User <: AbstractParameterNode allocated::Vector{Vector{Float64}} return_factor::Vector{Float64} min_level::Vector{Float64} - priorities::Vector{Int} record::@NamedTuple{ time::Vector{Float64}, subnetwork_id::Vector{Int}, @@ -507,7 +508,6 @@ struct User <: AbstractParameterNode allocated, return_factor, min_level, - priorities, record, ) else @@ -516,6 +516,19 @@ struct User <: AbstractParameterNode end end +""" +node_id: node ID of the AllocationTarget node +min_level: The minimum target level of the connected basin(s) +max_level: The maximum target level of the connected basin(s) +priority: If in a shortage state, the priority of the demand of the connected basin(s) +""" +struct AllocationTarget + node_id::Vector{NodeID} + min_level::Vector{LinearInterpolation} + max_level::Vector{LinearInterpolation} + priority::Vector{Int} +end + "Subgrid linearly interpolates basin levels." struct Subgrid basin_index::Vector{Int} @@ -558,6 +571,6 @@ struct Parameters{T, C1, C2} discrete_control::DiscreteControl pid_control::PidControl{T} user::User - lookup::Dict{Int, Symbol} + allocation_target::AllocationTarget subgrid::Subgrid end diff --git a/core/src/read.jl b/core/src/read.jl index 415cda76f..85abc3651 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -205,6 +205,10 @@ function initialize_allocation!(p::Parameters, config::Config)::Nothing (; allocation_network_ids, allocation_models, main_network_connections) = allocation allocation_network_ids_ = sort(collect(keys(graph[].node_ids))) + if isempty(allocation_network_ids_) + return nothing + end + errors = non_positive_allocation_network_id(graph) if errors error("Allocation network initialization failed.") @@ -646,8 +650,8 @@ function User(db::DB, config::Config)::User error("Problems encountered when parsing User static and time node IDs.") end - # All provided priorities - priorities = sort(unique(union(static.priority, time.priority))) + # All priorities used in the model + priorities = get_all_priorities(db, config) active = BitVector() min_level = Float64[] @@ -761,6 +765,31 @@ function User(db::DB, config::Config)::User ) end +function AllocationTarget(db::DB, config::Config)::AllocationTarget + static = load_structvector(db, config, AllocationTargetStaticV1) + time = load_structvector(db, config, AllocationTargetTimeV1) + + parsed_parameters, valid = parse_static_and_time( + db, + config, + "AllocationTarget"; + static, + time, + time_interpolatables = [:min_level, :max_level], + ) + + if !valid + error("Errors occurred when parsing AllocationTarget data.") + end + + return AllocationTarget( + NodeID.(NodeType.AllocationTarget, parsed_parameters.node_id), + parsed_parameters.min_level, + parsed_parameters.max_level, + parsed_parameters.priority, + ) +end + function Subgrid(db::DB, config::Config, basin::Basin)::Subgrid node_to_basin = Dict(node_id => index for (index, node_id) in enumerate(basin.node_id)) tables = load_structvector(db, config, BasinSubgridV1) @@ -815,6 +844,7 @@ function Parameters(db::DB, config::Config)::Parameters Int[], AllocationModel[], Vector{Tuple{NodeID, NodeID}}[], + get_all_priorities(db, config), Dict{Tuple{NodeID, NodeID}, Float64}(), Dict{Tuple{NodeID, NodeID}, Float64}(), (; @@ -845,6 +875,7 @@ function Parameters(db::DB, config::Config)::Parameters discrete_control = DiscreteControl(db, config) pid_control = PidControl(db, config, chunk_sizes) user = User(db, config) + allocation_target = AllocationTarget(db, config) basin = Basin(db, config, chunk_sizes) subgrid_level = Subgrid(db, config, basin) @@ -882,7 +913,7 @@ function Parameters(db::DB, config::Config)::Parameters discrete_control, pid_control, user, - Dict{Int, Symbol}(), + allocation_target, subgrid_level, ) diff --git a/core/src/schema.jl b/core/src/schema.jl index 5209152ce..0cc806824 100644 --- a/core/src/schema.jl +++ b/core/src/schema.jl @@ -23,6 +23,8 @@ @schema "ribasim.outlet.static" OutletStatic @schema "ribasim.user.static" UserStatic @schema "ribasim.user.time" UserTime +@schema "ribasim.allocationtarget.static" AllocationTargetStatic +@schema "ribasim.allocationtarget.time" AllocationTargetTime const delimiter = " / " tablename(sv::Type{SchemaVersion{T, N}}) where {T, N} = tablename(sv()) @@ -235,3 +237,18 @@ end min_level::Float64 priority::Int end + +@version AllocationTargetStaticV1 begin + node_id::Int + min_level::Float64 + max_level::Float64 + priority::Int +end + +@version AllocationTargetTimeV1 begin + node_id::Int + time::DateTime + min_level::Float64 + max_level::Float64 + priority::Int +end diff --git a/core/src/solve.jl b/core/src/solve.jl index d1383fc4c..4a8a04aa6 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -622,7 +622,7 @@ function formulate_du!( basin::Basin, storage::AbstractVector, )::Nothing - (; flow_vertical_dict, flow_vertical) = graph[] + (; flow_vertical) = graph[] flow_vertical = get_tmp(flow_vertical, storage) # loop over basins # subtract all outgoing flows diff --git a/core/src/util.jl b/core/src/util.jl index 6433b30a0..12c26d59f 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -620,22 +620,46 @@ function is_main_network(allocation_network_id::Int)::Bool return allocation_network_id == 1 end -function get_user_demand(user::User, node_id::NodeID, priority_idx::Int)::Float64 +function get_user_demand(p::Parameters, node_id::NodeID, priority_idx::Int)::Float64 + (; user, allocation) = p (; demand) = user user_idx = findsorted(user.node_id, node_id) - n_priorities = length(user.priorities) + n_priorities = length(allocation.priorities) return demand[(user_idx - 1) * n_priorities + priority_idx] end function set_user_demand!( - user::User, + p::Parameters, node_id::NodeID, priority_idx::Int, value::Float64, )::Nothing + (; user, allocation) = p (; demand) = user user_idx = findsorted(user.node_id, node_id) - n_priorities = length(user.priorities) + n_priorities = length(allocation.priorities) demand[(user_idx - 1) * n_priorities + priority_idx] = value return nothing end + +function get_all_priorities(db::DB, config::Config)::Vector{Int} + priorities = Set{Int}() + + # TODO: Is there a way to automatically grab all tables with a priority column? + for type in [UserStaticV1, UserTimeV1, AllocationTargetStaticV1, AllocationTargetTimeV1] + union!(priorities, load_structvector(db, config, type).priority) + end + return sort(unique(priorities)) +end + +function get_basin_priority(p::Parameters, node_id::NodeID)::Int + (; graph, allocation_target) = p + @assert node_id.type == NodeType.Basin + inneighbors_control = inneighbor_labels_type(graph, node_id, EdgeType.control) + if isempty(inneighbors_control) + return 0 + else + idx = findsorted(allocation_target.node_id, only(inneighbors_control)) + return allocation_target.priority[idx] + end +end diff --git a/core/src/validation.jl b/core/src/validation.jl index c84322fa1..caf0785c7 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -3,6 +3,7 @@ neighbortypes(nodetype::Symbol) = neighbortypes(Val(nodetype)) neighbortypes(::Val{:pump}) = Set((:basin, :fractional_flow, :terminal, :level_boundary)) neighbortypes(::Val{:outlet}) = Set((:basin, :fractional_flow, :terminal, :level_boundary)) neighbortypes(::Val{:user}) = Set((:basin, :fractional_flow, :terminal, :level_boundary)) +neighbortypes(::Val{:allocation_target}) = Set((:basin,)) neighbortypes(::Val{:basin}) = Set(( :linear_resistance, :tabulated_rating_curve, @@ -57,11 +58,12 @@ n_neighbor_bounds_flow(::Val{:Terminal}) = n_neighbor_bounds(1, typemax(Int), 0, n_neighbor_bounds_flow(::Val{:PidControl}) = n_neighbor_bounds(0, 0, 0, 0) n_neighbor_bounds_flow(::Val{:DiscreteControl}) = n_neighbor_bounds(0, 0, 0, 0) n_neighbor_bounds_flow(::Val{:User}) = n_neighbor_bounds(1, 1, 1, 1) +n_neighbor_bounds_flow(::Val{:AllocationTarget}) = n_neighbor_bounds(0, 0, 0, 0) n_neighbor_bounds_flow(nodetype) = error("'n_neighbor_bounds_flow' not defined for $nodetype.") n_neighbor_bounds_control(nodetype::Symbol) = n_neighbor_bounds_control(Val(nodetype)) -n_neighbor_bounds_control(::Val{:Basin}) = n_neighbor_bounds(0, 0, 0, typemax(Int)) +n_neighbor_bounds_control(::Val{:Basin}) = n_neighbor_bounds(0, 1, 0, typemax(Int)) n_neighbor_bounds_control(::Val{:LinearResistance}) = n_neighbor_bounds(0, 1, 0, 0) n_neighbor_bounds_control(::Val{:ManningResistance}) = n_neighbor_bounds(0, 1, 0, 0) n_neighbor_bounds_control(::Val{:TabulatedRatingCurve}) = n_neighbor_bounds(0, 1, 0, 0) @@ -75,6 +77,8 @@ n_neighbor_bounds_control(::Val{:PidControl}) = n_neighbor_bounds(0, 1, 1, 1) n_neighbor_bounds_control(::Val{:DiscreteControl}) = n_neighbor_bounds(0, 0, 1, typemax(Int)) n_neighbor_bounds_control(::Val{:User}) = n_neighbor_bounds(0, 0, 0, 0) +n_neighbor_bounds_control(::Val{:AllocationTarget}) = + n_neighbor_bounds(0, 0, 1, typemax(Int)) n_neighbor_bounds_control(nodetype) = error("'n_neighbor_bounds_control' not defined for $nodetype.") diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 1d2fbb4db..4a819dd16 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -1,5 +1,6 @@ @testitem "Allocation solve" begin using Ribasim: NodeID + using ComponentArrays: ComponentVector import SQLite import JuMP @@ -25,7 +26,8 @@ Ribasim.set_flow!(graph, NodeID(:FlowBoundary, 1), NodeID(:Basin, 2), 4.5) # Source flow allocation_model = p.allocation.allocation_models[1] - Ribasim.allocate!(p, allocation_model, 0.0) + u = ComponentVector(; storage = zeros(length(p.basin.node_id))) + Ribasim.allocate!(p, allocation_model, 0.0, u) F = allocation_model.problem[:F] @test JuMP.value(F[(NodeID(:Basin, 2), NodeID(:Basin, 6))]) ≈ 0.0 @@ -42,9 +44,9 @@ # Test getting and setting user demands (; user) = p - Ribasim.set_user_demand!(user, NodeID(:User, 11), 2, Float64(π)) + Ribasim.set_user_demand!(p, NodeID(:User, 11), 2, Float64(π)) @test user.demand[4] ≈ π - @test Ribasim.get_user_demand(user, NodeID(:User, 11), 2) ≈ π + @test Ribasim.get_user_demand(p, NodeID(:User, 11), 2) ≈ π end @testitem "Allocation objective: quadratic absolute" skip = true begin @@ -118,12 +120,12 @@ end problem = model.integrator.p.allocation.allocation_models[1].problem objective = JuMP.objective_function(problem) @test objective isa JuMP.AffExpr # Affine expression - @test :F_abs in keys(problem.obj_dict) + @test :F_abs_user in keys(problem.obj_dict) F = problem[:F] - F_abs = problem[:F_abs] + F_abs_user = problem[:F_abs_user] - @test objective.terms[F_abs[NodeID(:User, 5)]] == 1.0 - @test objective.terms[F_abs[NodeID(:User, 6)]] == 1.0 + @test objective.terms[F_abs_user[NodeID(:User, 5)]] == 1.0 + @test objective.terms[F_abs_user[NodeID(:User, 6)]] == 1.0 end @testitem "Allocation objective: linear relative" begin @@ -142,12 +144,12 @@ end problem = model.integrator.p.allocation.allocation_models[1].problem objective = JuMP.objective_function(problem) @test objective isa JuMP.AffExpr # Affine expression - @test :F_abs in keys(problem.obj_dict) + @test :F_abs_user in keys(problem.obj_dict) F = problem[:F] - F_abs = problem[:F_abs] + F_abs_user = problem[:F_abs_user] - @test objective.terms[F_abs[NodeID(:User, 5)]] == 1.0 - @test objective.terms[F_abs[NodeID(:User, 6)]] == 1.0 + @test objective.terms[F_abs_user[NodeID(:User, 5)]] == 1.0 + @test objective.terms[F_abs_user[NodeID(:User, 6)]] == 1.0 end @testitem "Allocation with controlled fractional flow" begin @@ -243,9 +245,9 @@ end # support absolute value expressions in the objective function allocation_model_main_network = Ribasim.get_allocation_model(p, 1) problem = allocation_model_main_network.problem - @test problem[:F_abs].axes[1] == NodeID.(:Pump, [11, 24, 38]) - @test problem[:abs_positive].axes[1] == NodeID.(:Pump, [11, 24, 38]) - @test problem[:abs_negative].axes[1] == NodeID.(:Pump, [11, 24, 38]) + @test problem[:F_abs_user].axes[1] == NodeID.(:Pump, [11, 24, 38]) + @test problem[:abs_positive_user].axes[1] == NodeID.(:Pump, [11, 24, 38]) + @test problem[:abs_negative_user].axes[1] == NodeID.(:Pump, [11, 24, 38]) # In each subnetwork, the connection from the main network to the subnetwork is # interpreted as a source @@ -260,6 +262,7 @@ end @testitem "allocation with main network optimization problem" begin using SQLite using Ribasim: NodeID + using ComponentArrays: ComponentVector using JuMP toml_path = normpath( @@ -273,13 +276,14 @@ end p = Ribasim.Parameters(db, cfg) close(db) - (; allocation, user, graph) = p + (; allocation, user, graph, basin) = p (; allocation_models, subnetwork_demands, subnetwork_allocateds) = allocation t = 0.0 # Collecting demands + u = ComponentVector(; storage = zeros(length(basin.node_id))) for allocation_model in allocation_models[2:end] - Ribasim.allocate!(p, allocation_model, t; collect_demands = true) + Ribasim.allocate!(p, allocation_model, t, u; collect_demands = true) end @test subnetwork_demands[(NodeID(:Basin, 2), NodeID(:Pump, 11))] ≈ [4.0, 4.0, 0.0] @@ -291,19 +295,20 @@ end # containing subnetworks as users allocation_model = allocation_models[1] (; problem) = allocation_model - Ribasim.allocate!(p, allocation_model, t) + Ribasim.allocate!(p, allocation_model, t, u) # Main network objective function objective = JuMP.objective_function(problem) objective_variables = keys(objective.terms) - F_abs = problem[:F_abs] - @test F_abs[NodeID(:Pump, 11)] ∈ objective_variables - @test F_abs[NodeID(:Pump, 24)] ∈ objective_variables - @test F_abs[NodeID(:Pump, 38)] ∈ objective_variables + F_abs_user = problem[:F_abs_user] + @test F_abs_user[NodeID(:Pump, 11)] ∈ objective_variables + @test F_abs_user[NodeID(:Pump, 24)] ∈ objective_variables + @test F_abs_user[NodeID(:Pump, 38)] ∈ objective_variables # Running full allocation algorithm Ribasim.set_flow!(graph, NodeID(:FlowBoundary, 1), NodeID(:Basin, 2), 4.5) - Ribasim.update_allocation!((; p, t)) + u = ComponentVector(; storage = zeros(length(p.basin.node_id))) + Ribasim.update_allocation!((; p, t, u)) @test subnetwork_allocateds[NodeID(:Basin, 2), NodeID(:Pump, 11)] ≈ [4.0, 0.49500000, 0.0] @@ -314,3 +319,54 @@ end @test user.allocated[2] ≈ [4.0, 0.0, 0.0] @test user.allocated[7] ≈ [0.001, 0.0, 0.0] end + +@testitem "Allocation level control" begin + toml_path = + normpath(@__DIR__, "../../generated_testmodels/allocation_target/ribasim.toml") + @test ispath(toml_path) + model = Ribasim.run(toml_path) + + storage = Ribasim.get_storages_and_levels(model).storage[1, :] + t = Ribasim.timesteps(model) + + p = model.integrator.p + (; user, graph, allocation, basin, allocation_target) = p + + d = user.demand_itp[1][2](0) + ϕ = 1e-3 # precipitation + q = Ribasim.get_flow( + graph, + Ribasim.NodeID(Ribasim.NodeType.FlowBoundary, 1), + Ribasim.NodeID(Ribasim.NodeType.Basin, 2), + 0, + ) + A = basin.area[1][1] + l_max = allocation_target.max_level[1](0) + Δt_allocation = allocation.allocation_models[1].Δt_allocation + + # Until the first allocation solve, the user abstracts fully + pre_allocation = t .<= Δt_allocation + u_pre_allocation(τ) = storage[1] + (q + ϕ - d) * τ + @test storage[pre_allocation] ≈ u_pre_allocation.(t[pre_allocation]) rtol = 1e-4 + + # Until the basin is at its maximum level, the user does not abstract + basin_filling = @. ~pre_allocation && (storage <= A * l_max) + fill_start_idx = findlast(pre_allocation) + u_filling(τ) = storage[fill_start_idx] + (q + ϕ) * (τ - t[fill_start_idx]) + @test storage[basin_filling] ≈ u_filling.(t[basin_filling]) rtol = 1e-4 + + # After the basin has reached its maximum level, the user abstracts fully again + precipitation = eachindex(storage) .<= argmax(storage) + after_filling = @. ~pre_allocation && ~basin_filling && precipitation + fill_stop_idx = findfirst(after_filling) + u_after_filling(τ) = storage[fill_stop_idx] + (q + ϕ - d) * (τ - t[fill_stop_idx]) + @test storage[after_filling] ≈ u_after_filling.(t[after_filling]) rtol = 1e-4 + + # After precipitation stops, the user still abstracts from the basin so the storage decreases + storage_reduction = @. ~precipitation && t <= 1.8e6 + storage_reduction_start = findfirst(storage_reduction) + u_storage_reduction(τ) = + storage[storage_reduction_start] + (q - d) * (τ - t[storage_reduction_start]) + @test storage[storage_reduction] ≈ u_storage_reduction.(t[storage_reduction]) rtol = + 1e-4 +end diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index a97773330..b9d821348 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -54,11 +54,19 @@ $$ \phi_i(t), \quad \forall i \in B_S. $$ -Secondly, there is the available water in each basin above the minimum level $l_{\min,i}$ corresponding to a minimum storage $s_{\min,i}$: +We consider fluxes into the basin to be positive and out of the basin to be negative. + +Secondly, there is either a supply or demand from the storage in the basin. Given a minimum level $\ell_{\min, i}$ and a maximum level $\ell_{\max, i}$ which correspond to a minimum storage $s_{\min, i}$ and maximum storage $s_{\max, i}$ respectively, we get a flow supply of +$$ + F^{\text{basin out}}_{\max, i} = \max\left(0.0, \frac{u_i(t)-s_{\max,i}}{\Delta t_{\text{alloc}}} + \phi_i(t)\right) $$ - u_i(t)-s_{\min,i}, \quad \forall i \in B_S. + +and a demand of $$ -Note that this value can be negative, which we interpret as a demand from the basin. + d^p_i = \max\left(0.0, \frac{s_{\min,i} - u_i(t)}{\Delta t_{\text{alloc}}} - \phi_i(t)\right), +$$ + +for all $i \in B_S$. Here $\Delta t_{\text{alloc}}$ is the simulated time between two consecutive allocation solves. Note that the basin demand has only a single priority, so for other priorities this demand is $0$. ### Constraining factors @@ -92,49 +100,50 @@ for the set of in-neighbors and out-neighbors of a node in the allocation networ Each edge in the allocation network has an associated capacity. These capacities are collected in the sparse capacity matrix $C_S \in \overline{\mathbb{R}}_{\ge 0}^{n'\times n'}$ where $n' = \#V'_S$ is the number of nodes in the allocation network. The capacities can be infinite, if there is nothing in the model constraining the capacity of the edge. -The capacities are determined in 4 different ways: +The capacities are determined in different ways: - If an edge does not exist in the allocation network, i.e. $(i,j) \notin E_S$ for certain $1 \le i,j\le n'$, then $(C_S)_{i,j} = 0$; - The capacity of the edge $e \in E_S$ is given by the smallest `max_flow_rate` of the nodes along the equivalent edges in the subnetwork. If there are no nodes with a `max_flow_rate`, the edge capacity is infinite; - If the edge is a source, the capacity of the edge is given by the flow rate of that source. +There are also capacities $C^B_S \in \mathbb{R}^b$ where $b = \# B_S$ is the number of basins, for the flow supplied by basins. + # The optimization problem The optimization problem for a subnetwork is a linear optimization problem consisting of an objective function with associated constraints on a set of variables, all of which are introduced below. ## The optimization variables -There are 2 types of variable whose value has to be determined to solve the allocation problem: +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 allocations to the basins -$$ - A^\text{basin}_{i} \ge 0, \quad B_S. -$$ - -:::{.callout-note} -Currently the basin allocations are not taken into account in the implementation. -::: +- 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 optimization objective -The goal of allocation is to get the flow to the users as close as possible to their demand. To achieve this, the following objectives are supported: +The goal of allocation is to get the flow to nodes with demands as close as possible to these demands. To achieve this, a sum error of terms is minimized. The form of these error terms is determined by the choice of one of the supported objective function types, which are discussed further below. + +$$ + \min E_{\text{user}} + E_{\text{basin}} +$$ + +### User demands - `quadratic_absolute`: $$ - \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left( F_{ij} - d_j^p(t)\right)^2 + E_{\text{user}} = \sum_{(i,j)\in E_S\;:\; i\in U_S} \left( F_{ij} - d_j^p(t)\right)^2 $$ - `quadratic_relative`: $$ - \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left( 1 - \frac{F_{ij}}{d_j^p(t)}\right)^2 + E_{\text{user}} = \sum_{(i,j)\in E_S\;:\; i\in U_S} \left( 1 - \frac{F_{ij}}{d_j^p(t)}\right)^2 $$ - `linear_absolute` (default): $$ - \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left| F_{ij} - d_j^p(t)\right| + E_{\text{user}} = \sum_{(i,j)\in E_S\;:\; i\in U_S} \left| F_{ij} - d_j^p(t)\right| $$ - `linear_relative`: $$ - \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left|1 - \frac{F_{ij}}{d_j^p(t)}\right| + E_{\text{user}} = \sum_{(i,j)\in E_S\;:\; i\in U_S} \left|1 - \frac{F_{ij}}{d_j^p(t)}\right| $$ :::{.callout-note} @@ -156,16 +165,45 @@ The absolute value applied here is not supported in a linear programming context In the future new optimization objectives will be introduced, for demands of basins and priorities over sources. These will be used in combination with the above, in the form of goal programming. ::: +### Basin demands + +- `quadratic_absolute`: +$$ + E_{\text{basin}} = \sum_{i \in B_S} \left( F_i^\text{basin in} - d_i^p(t)\right)^2 +$$ +- `quadratic_relative`: +$$ + E_{\text{basin}} = \sum_{i \in B_S} \left( 1 - \frac{F_i^\text{basin in}}{d_i^p(t)}\right)^2 +$$ +- `linear_absolute` (default): +$$ + E_{\text{basin}} = \sum_{i \in B_S} \left| F_i^\text{basin in} - d_i^p(t)\right| +$$ +- `linear_relative`: +$$ + E_{\text{basin}} = \sum_{i \in B_S} \left|1 - \frac{F_i^\text{basin in}}{d_i^p(t)}\right| +$$ + ## The optimization constraints - Flow conservation: For the basins in the allocation network we have that $$ - \sum_{j=1}^{n'} F_{kj} = \sum_{i=1}^{n'} F_{ik}, \quad \forall k \in B_S. + F^\text{basin in}_k + \sum_{j=1}^{n'} F_{kj} = F^\text{basin out}_k + \sum_{i=1}^{n'} F_{ik}, \quad \forall k \in B_S . $$ {#eq-flowconservationconstraint} +Note that we do not require equality here; in the allocation we do not mind that excess flow is 'forgotten' if it cannot contribute to the allocation to the users. + +:::{.callout-note} +In @eq-flowconservationconstraint, the placement of the basin flows might seem counter-intuitive. Think of the basin storage as a separate node connected to the basin node. +::: + - Capacity: the flows over the edges are positive and bounded by the edge capacity: $$ F_{ij} \le \left(C_S\right)_{ij}, \quad \forall(i,j) \in E_S. $$ {#eq-capacityconstraint} -By the definition of $C_S$ this also includes the source flows. +By the definition of $C_S$ this also includes the source flows. The same holds for the basin outflows: + +$$ + F^{\text{basin out}}_{i} \le F^{\text{basin out}}_{\max, i}, \quad \forall i \in B_S. +$$ :::{.callout-note} When performing subnetwork demand collection, these capacities are set to $\infty$ for edges which connect the main network to a subnetwork. @@ -232,9 +270,11 @@ The following is an example of an optimization problem for the example shown [he using Ribasim using Ribasim: NodeID using SQLite +using ComponentArrays: ComponentVector toml_path = normpath(@__DIR__, "../../generated_testmodels/allocation_example/ribasim.toml") p = Ribasim.Model(toml_path).integrator.p +u = ComponentVector(; storage = zeros(length(p.basin.node_id))) allocation_model = p.allocation.allocation_models[1] t = 0.0 @@ -244,7 +284,7 @@ Ribasim.set_flow!(p.graph, NodeID(:FlowBoundary, 1), NodeID(:Basin, 2), 1.0) Ribasim.adjust_source_capacities!(allocation_model, p, priority_idx) Ribasim.adjust_edge_capacities!(allocation_model, p, priority_idx) -Ribasim.set_objective_priority!(allocation_model, p, t, priority_idx) +Ribasim.set_objective_priority!(allocation_model, p, u, t, priority_idx) println(p.allocation.allocation_models[1].problem) ``` diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index 1cc600b8e..28cbcca9c 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -198,6 +198,9 @@ name it must have in the database if it is stored there. - User: sets water usage demands at certain priorities - `User / static`: demands - `User / time`: dynamic demands +- AllocationTarget: Indicates minimum and maximum target level of connected basins for allocation + - `AllocationTarget / static`: static target levels + - `AllocationTarget / time`: dynamic target levels - Terminal: Water sink without state or properties - `Terminal / static`: - (only node IDs) - DiscreteControl: Set parameters of other nodes based on model state conditions (e.g. basin level) @@ -342,37 +345,6 @@ Water levels beyond the last `basin_level` are linearly extrapolated. Note that the interpolation to subgrid water level is not constrained by any water balance within Ribasim. Generally, to create physically meaningful subgrid water levels, the subgrid table must be parametrized properly such that the spatially integrated water volume of the subgrid elements agrees with the total storage volume of the basin. -## Basin results - -The basin table contains results of the storage and level of each basin at every solver -timestep. The initial condition is also written to the file. - -column | type | unit --------- | -------- | ---- -time | DateTime | - -node_id | Int | - -storage | Float64 | $m^3$ -level | Float64 | $m$ - -The table is sorted by time, and per time it is sorted by `node_id`. - -## Flow results - -The flow table contains results of the flow on every edge in the model, for each solver -timestep. - -column | type | unit -------------- | ------------------- | ---- -time | DateTime | - -edge_id | Union{Int, Missing} | - -from_node_id | Int | - -to_node_id | Int | - -flow | Float64 | $m^3 s^{-1}$ - -The table is sorted by time, and per time the same `edge_id` order is used, though not sorted. -Flows that are added to the model at a node, have a missing `edge_id`, and identical `from_node_id` and `to_node_id`. -Flows out of the model always have a negative sign, and additions a positive sign. - # FractionalFlow {#sec-fractional-flow} Lets a fraction (in [0,1]) of the incoming flow trough. @@ -470,8 +442,8 @@ node_id | Int | - | sorted active | Bool | - | (optional, default true) demand | Float64 | $m^3 s^{-1}$ | - return_factor | Float64 | - | between [0 - 1] -min_level | Float64 | $m$ | (optional) -priority | Int | - | sorted per node id +min_level | Float64 | $m$ | - +priority | Int | - | positive, sorted per node id ## User / time @@ -486,29 +458,38 @@ Note that a `node_id` can be either in this table or in the static one, but not column | type | unit | restriction ------------- | -------- | ------------ | ----------- node_id | Int | - | sorted -priority | Int | - | sorted per node id +priority | Int | - | positive, sorted per node id time | DateTime | - | sorted per priority per node id demand | Float64 | $m^3 s^{-1}$ | - return_factor | Float64 | - | between [0 - 1] -min_level | Float64 | $m$ | (optional) +min_level | Float64 | $m$ | - -## Allocation results +# AllocationTarget {#sec-allocation_target} -The allocation table contains a record of allocation results: when it happened, for which user node, in which allocation network, and what the demand, allocated flow and abstracted flow were. +An `AllocationTarget` node associates a minimum and a maximum level with connected basins to be used by the allocation algorithm. +Below the minimum level the basin has a demand of the supplied priority, +above the maximum level the basin acts as a source, used by all nodes with demands in order of priority. +The same `AllocationTarget` node can be used for basins in different subnetworks. -column | type ---------------| ----- -time | DateTime -subnetwork_id | Int -user_node_id | Int -priority | Int -demand | Float64 -allocated | Float64 -abstracted | Float64 +column | type | unit | restriction +------------- | ------- | ------------ | ----------- +node_id | Int | - | sorted +min_level | Float64 | $m$ | - +max_level | Float64 | $m$ | - +priority | Int | - | positive -::: {.callout-note} -Currently the stored demand and abstraction rate are those at the allocation timepoint (and the abstraction rate is based on the previous allocation optimization). In the future these will be an average over the previous allocation timestep. -::: +## AllocationTarget / time + +This table is the transient form of the `AllocationTarget` table, in which time-dependent minimum and maximum levels can be supplied. +Similar to the static version, only a single priority per `AllocationTarget` node can be provided. + +column | type | unit | restriction +------------- | ------- | ------------ | ----------- +node_id | Int | - | sorted +time | DateTime | - | sorted per node id +min_level | Float64 | $m$ | - +max_level | Float64 | $m$ | - +priority | Int | - | positive # LevelBoundary {#sec-level-boundary} @@ -642,17 +623,6 @@ node_id | Int | - | sorted control_state | String | - | - truth_state | String | - | Consists of the characters "T" (true), "F" (false), "U" (upcrossing), "D" (downcrossing) and "*" (any), sorted per node_id -## DiscreteControl results - -The control table contains a record of each change of control state: when it happened, which control node was involved, to which control state it changed and based on which truth state. - -column | type ---------------- | ---------- -time | DateTime -control_node_id | Int -truth_state | String -control_state | String - # PidControl {#sec-pid-control} The PidControl node controls the level in a basin by continuously controlling the flow rate of a connected pump or outlet. See also [PID controller](https://en.wikipedia.org/wiki/PID_controller). When A PidControl node is made inactive, the node under its control retains the last flow rate value, and the error integral is reset to 0. @@ -689,3 +659,82 @@ target | Float64 | $m$ | - proportional | Float64 | $s^{-1}$ | - integral | Float64 | $s^{-2}$ | - derivative | Float64 | - | - + +# Results + +## Basin - `basin.arrow` + +The basin table contains results of the storage and level of each basin at every solver +timestep. The initial condition is also written to the file. + +column | type | unit +-------- | -------- | ---- +time | DateTime | - +node_id | Int | - +storage | Float64 | $m^3$ +level | Float64 | $m$ + +The table is sorted by time, and per time it is sorted by `node_id`. + +## Flow - `flow.arrow` + +The flow table contains results of the flow on every edge in the model, for each solver +timestep. + +column | type | unit +------------- | ------------------- | ---- +time | DateTime | - +edge_id | Union{Int, Missing} | - +from_node_id | Int | - +to_node_id | Int | - +flow | Float64 | $m^3 s^{-1}$ + +The table is sorted by time, and per time the same `edge_id` order is used, though not sorted. +Flows that are added to the model at a node, have a missing `edge_id`, and identical `from_node_id` and `to_node_id`. +Flows out of the model always have a negative sign, and additions a positive sign. + +## DiscreteControl - `control.arrow` + +The control table contains a record of each change of control state: when it happened, which control node was involved, to which control state it changed and based on which truth state. + +column | type +--------------- | ---------- +time | DateTime +control_node_id | Int +truth_state | String +control_state | String + +## Allocation - `allocation.arrow` + +The allocation table contains a record of allocation results: when it happened, for which user node, in which allocation network, and what the demand, allocated flow and abstracted flow were. + +column | type +--------------| ----- +time | DateTime +subnetwork_id | Int +user_node_id | Int +priority | Int +demand | Float64 +allocated | Float64 +abstracted | Float64 + +::: {.callout-note} +Currently the stored demand and abstraction rate are those at the allocation timepoint (and the abstraction rate is based on the previous allocation optimization). In the future these will be an average over the previous allocation timestep. +::: + +## Allocation flow - `allocation_flow.arrow` + +The allocation flow table contains results of the optimized allocation flow on every edge in the model that is part of a subnetwork, for each time an optimization problem is solved (see also [here](allocation.qmd#the-high-level-algorithm)). +If in the model a main network and subnetwork(s) are specified, there are 2 different +types of optimization for the subnetwork: collecting its total demand per priority (for allocating flow from the main network to the subnetwork), and allocating flow within the subnetwork. The column `collect_demands` provides the distinction between these two optimization types. + +column | type +----------------| ----- +time | DateTime +edge_id | Int +from_node_id | Int +to_node_id | Int +subnetwork_id | Int +priority | Int +flow | Float64 +collect_demands | Bool diff --git a/python/ribasim/ribasim/__init__.py b/python/ribasim/ribasim/__init__.py index c788eee55..85322bdce 100644 --- a/python/ribasim/ribasim/__init__.py +++ b/python/ribasim/ribasim/__init__.py @@ -4,6 +4,7 @@ from ribasim import utils from ribasim.config import ( Allocation, + AllocationTarget, Basin, Compression, DiscreteControl, @@ -29,6 +30,7 @@ __all__ = [ "Allocation", + "AllocationTarget", "Basin", "DiscreteControl", "Compression", diff --git a/python/ribasim/ribasim/config.py b/python/ribasim/ribasim/config.py index 2189366e9..8c046e1f8 100644 --- a/python/ribasim/ribasim/config.py +++ b/python/ribasim/ribasim/config.py @@ -7,6 +7,8 @@ # These schemas are autogenerated from ribasim.schemas import ( + AllocationTargetStaticSchema, + AllocationTargetTimeSchema, BasinProfileSchema, BasinStateSchema, BasinStaticSchema, @@ -136,6 +138,17 @@ class User(NodeModel): ) +class AllocationTarget(NodeModel): + static: TableModel[AllocationTargetStaticSchema] = Field( + default_factory=TableModel[AllocationTargetStaticSchema], + json_schema_extra={"sort_keys": ["node_id", "priority"]}, + ) + time: TableModel[AllocationTargetTimeSchema] = Field( + default_factory=TableModel[AllocationTargetTimeSchema], + json_schema_extra={"sort_keys": ["node_id", "priority", "time"]}, + ) + + class FlowBoundary(NodeModel): static: TableModel[FlowBoundaryStaticSchema] = Field( default_factory=TableModel[FlowBoundaryStaticSchema], diff --git a/python/ribasim/ribasim/geometry/node.py b/python/ribasim/ribasim/geometry/node.py index 3272f132f..cb6044229 100644 --- a/python/ribasim/ribasim/geometry/node.py +++ b/python/ribasim/ribasim/geometry/node.py @@ -205,6 +205,7 @@ def plot(self, ax=None, zorder=None) -> Any: "DiscreteControl": "*", "PidControl": "x", "User": "s", + "AllocationTarget": "o", "": "o", } @@ -222,6 +223,7 @@ def plot(self, ax=None, zorder=None) -> Any: "DiscreteControl": "k", "PidControl": "k", "User": "g", + "AllocationTarget": "k", "": "k", } assert self.df is not None diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py index f485477a3..dfc1541fd 100644 --- a/python/ribasim/ribasim/model.py +++ b/python/ribasim/ribasim/model.py @@ -18,6 +18,7 @@ from ribasim.config import ( Allocation, + AllocationTarget, Basin, DiscreteControl, FlowBoundary, @@ -168,6 +169,7 @@ class Model(FileModel): logging: Logging = Logging() allocation: Allocation = Field(default_factory=Allocation) + allocation_target: AllocationTarget = Field(default_factory=AllocationTarget) basin: Basin = Field(default_factory=Basin) fractional_flow: FractionalFlow = Field(default_factory=FractionalFlow) level_boundary: LevelBoundary = Field(default_factory=LevelBoundary) diff --git a/python/ribasim/ribasim/schemas.py b/python/ribasim/ribasim/schemas.py index 890649b07..44e9dacd0 100644 --- a/python/ribasim/ribasim/schemas.py +++ b/python/ribasim/ribasim/schemas.py @@ -11,6 +11,21 @@ class Config: coerce = True +class AllocationTargetStaticSchema(_BaseSchema): + node_id: Series[int] = pa.Field(nullable=False) + min_level: Series[float] = pa.Field(nullable=False) + max_level: Series[float] = pa.Field(nullable=False) + priority: Series[int] = pa.Field(nullable=False) + + +class AllocationTargetTimeSchema(_BaseSchema): + node_id: Series[int] = pa.Field(nullable=False) + time: Series[Timestamp] = pa.Field(nullable=False) + min_level: Series[float] = pa.Field(nullable=False) + max_level: Series[float] = pa.Field(nullable=False) + priority: Series[int] = pa.Field(nullable=False) + + class BasinProfileSchema(_BaseSchema): node_id: Series[int] = pa.Field(nullable=False, default=0) area: Series[float] = pa.Field(nullable=False) diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index 4616d2746..107d6c3a5 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -7,6 +7,7 @@ import ribasim_testmodels from ribasim_testmodels.allocation import ( allocation_example_model, + allocation_target_model, fractional_flow_subnetwork_model, looped_subnetwork_model, main_network_with_subnetworks_model, @@ -54,6 +55,7 @@ __all__ = [ "allocation_example_model", + "allocation_target_model", "backwater_model", "basic_arrow_model", "basic_model", diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 9d7cb4fe5..285fa46e0 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -1015,7 +1015,9 @@ def allocation_example_model(): def main_network_with_subnetworks_model(): - # Set ip the nodes: + """Generate a model which consists of a main network and multiple connected subnetworks.""" + + # Set up the nodes: xy = np.array( [ (0.0, -1.0), @@ -1638,3 +1640,122 @@ def main_network_with_subnetworks_model(): ) return model + + +def allocation_target_model(): + # Set up the nodes: + xy = np.array( + [ + (0.0, 0.0), # 1: FlowBoundary + (1.0, 0.0), # 2: Basin + (2.0, 0.0), # 3: User + (1.0, -1.0), # 4: AllocationTarget + (2.0, -1.0), # 5: Basin + ] + ) + node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1]) + + node_type = ["FlowBoundary", "Basin", "User", "AllocationTarget", "Basin"] + + # Make sure the feature id starts at 1: explicitly give an index. + node = ribasim.Node( + df=gpd.GeoDataFrame( + data={ + "type": node_type, + "subnetwork_id": 5 * [2], + }, + index=pd.Index(np.arange(len(xy)) + 1, name="fid"), + geometry=node_xy, + crs="EPSG:28992", + ) + ) + + # Setup the edges: + from_id = np.array([1, 2, 4, 3, 4]) + to_id = np.array([2, 3, 2, 5, 5]) + edge_type = ["flow", "flow", "control", "flow", "control"] + allocation_network_id = [1, None, None, None, None] + + lines = node.geometry_from_connectivity(from_id.tolist(), to_id.tolist()) + edge = ribasim.Edge( + df=gpd.GeoDataFrame( + data={ + "from_node_id": from_id, + "to_node_id": to_id, + "edge_type": edge_type, + "allocation_network_id": allocation_network_id, + }, + geometry=lines, + crs="EPSG:28992", + ) + ) + + # Setup basin + profile = pd.DataFrame( + data={"node_id": [2, 2, 5, 5], "area": 1e3, "level": [0.0, 1.0, 0.0, 1.0]} + ) + static = pd.DataFrame( + data={ + "node_id": [5], + "drainage": 0.0, + "potential_evaporation": 0.0, + "infiltration": 0.0, + "precipitation": 0.0, + "urban_runoff": 0.0, + } + ) + time = pd.DataFrame( + data={ + "node_id": 2, + "time": ["2020-01-01 00:00:00", "2020-01-16 00:00:00"], + "drainage": 0.0, + "potential_evaporation": 0.0, + "infiltration": 0.0, + "precipitation": [1e-6, 0.0], + "urban_runoff": 0.0, + }, + ) + + state = pd.DataFrame(data={"node_id": [2, 5], "level": 0.5}) + basin = ribasim.Basin(profile=profile, static=static, time=time, state=state) + + # Setup flow boundary + flow_boundary = ribasim.FlowBoundary( + static=pd.DataFrame(data={"node_id": [1], "flow_rate": 1e-3}) + ) + + # Setup allocation level control + allocation_target = ribasim.AllocationTarget( + static=pd.DataFrame( + data={"node_id": [4], "priority": 1, "min_level": 1.0, "max_level": 1.5} + ) + ) + + # Setup user + user = ribasim.User( + static=pd.DataFrame( + data={ + "node_id": [3], + "priority": [2], + "demand": [1.5e-3], + "return_factor": [0.2], + "min_level": [0.2], + } + ) + ) + + # Setup allocation + allocation = ribasim.Allocation(use_allocation=True, timestep=1e5) + + model = ribasim.Model( + network=ribasim.Network(node=node, edge=edge), + basin=basin, + flow_boundary=flow_boundary, + allocation_target=allocation_target, + user=user, + allocation=allocation, + starttime="2020-01-01 00:00:00", + endtime="2020-02-01 00:00:00", + ) + + return model diff --git a/ribasim_qgis/core/nodes.py b/ribasim_qgis/core/nodes.py index 3cfb49be4..e8b77a39a 100644 --- a/ribasim_qgis/core/nodes.py +++ b/ribasim_qgis/core/nodes.py @@ -215,6 +215,11 @@ def renderer(self) -> QgsCategorizedSymbolRenderer: "DiscreteControl": (QColor("black"), "DiscreteControl", shape.Star), "PidControl": (QColor("black"), "PidControl", shape.Cross2), "User": (QColor("green"), "User", shape.Square), + "AllocationTarget": ( + QColor("black"), + "AllocationTarget", + shape.Circle, + ), # All other nodes, or incomplete input "": (QColor("white"), "", shape.Circle), } @@ -780,6 +785,43 @@ def attributes(cls) -> list[QgsField]: ] +class AllocationTargetStatic(Input): + @classmethod + def input_type(cls) -> str: + return "AllocationTarget / static" + + @classmethod + def geometry_type(cls) -> str: + return "No Geometry" + + @classmethod + def attributes(cls) -> list[QgsField]: + return [ + QgsField("node_id", QVariant.Int), + QgsField("target_level", QVariant.Double), + QgsField("priority", QVariant.Int), + ] + + +class AllocationTargetTime(Input): + @classmethod + def input_type(cls) -> str: + return "AllocationTarget / time" + + @classmethod + def geometry_type(cls) -> str: + return "No Geometry" + + @classmethod + def attributes(cls) -> list[QgsField]: + return [ + QgsField("node_id", QVariant.Int), + QgsField("time", QVariant.DateTime), + QgsField("target_level", QVariant.Double), + QgsField("priority", QVariant.Int), + ] + + NODES: dict[str, type[Input]] = { cls.input_type(): cls # type: ignore[type-abstract] # mypy doesn't see that all classes are concrete. for cls in Input.__subclasses__()