From 6910f4372fb2b300db65d09dbf2f31ff93a6dbe3 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 24 Apr 2024 10:15:14 +0200 Subject: [PATCH 1/4] Use mean source flows in allocation --- core/src/allocation_optim.jl | 7 ++++--- core/src/callback.jl | 25 +++++++++++++++++++++++-- core/src/graph.jl | 26 ++++++++++++++++---------- core/src/parameter.jl | 13 ++++++++----- core/src/read.jl | 23 ++++++++++++++++------- 5 files changed, 67 insertions(+), 27 deletions(-) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 75dc66c12..59ae9dc7c 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -254,7 +254,8 @@ function set_initial_capacities_source!( p::Parameters, )::Nothing (; problem) = allocation_model - (; graph) = p + (; graph, allocation) = p + (; mean_source_flows) = allocation (; subnetwork_id) = allocation_model source_constraints = problem[:source] main_network_source_edges = get_main_network_connections(p, subnetwork_id) @@ -264,8 +265,8 @@ function set_initial_capacities_source!( if graph[edge...].subnetwork_id_source == subnetwork_id # If it is a source edge for this allocation problem if edge ∉ main_network_source_edges - # Reset the source to the current flow from the physical layer. - source_capacity = get_flow(graph, edge..., 0) + # Reset the source to the averaged flow over the last allocation period + source_capacity = mean_source_flows[edge][] JuMP.set_normalized_rhs( source_constraints[edge], # It is assumed that the allocation procedure does not have to be differentiated. diff --git a/core/src/callback.jl b/core/src/callback.jl index 51dd3d029..b22b8901e 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -91,7 +91,7 @@ Integrate flows over the last timestep """ function integrate_flows!(u, t, integrator)::Nothing (; p, dt) = integrator - (; graph, user_demand, basin) = p + (; graph, user_demand, basin, allocation) = p (; flow, flow_dict, flow_prev, flow_integrated) = graph[] (; vertical_flux, vertical_flux_prev, vertical_flux_integrated, vertical_flux_bmi) = basin @@ -106,12 +106,21 @@ function integrate_flows!(u, t, integrator)::Nothing @. vertical_flux_integrated += 0.5 * (vertical_flux + vertical_flux_prev) * dt @. vertical_flux_bmi += 0.5 * (vertical_flux + vertical_flux_prev) * dt + # UserDemand realized flows for BMI for (i, id) in enumerate(user_demand.node_id) src_id = inflow_id(graph, id) flow_idx = flow_dict[src_id, id] user_demand.realized_bmi[i] += 0.5 * (flow[flow_idx] + flow_prev[flow_idx]) * dt end + # Allocation source flows + for (edge, value) in allocation.mean_source_flows + value[] += + 0.5 * + (get_flow(graph, edge..., 0) + get_flow(graph, edge..., 0; prev = true)) * + dt + end + copyto!(flow_prev, flow) copyto!(vertical_flux_prev, vertical_flux) return nothing @@ -446,7 +455,14 @@ end function update_allocation!(integrator)::Nothing (; p, t, u) = integrator (; allocation) = p - (; allocation_models) = allocation + (; allocation_models, mean_source_flows) = allocation + (; Δt_allocation) = allocation_models[1] + + # Divide by the allocation Δt to obtain the mean flows + # from the integrated flows + for value in values(mean_source_flows) + value[] /= Δt_allocation + end # If a main network is present, collect demands of subnetworks if has_main_network(allocation) @@ -462,6 +478,11 @@ function update_allocation!(integrator)::Nothing for allocation_model in allocation_models allocate!(p, allocation_model, t, u, OptimizationType.allocate) end + + # Reset the mean source flows + for value in values(mean_source_flows) + value[] = 0.0 + end end "Load updates from 'TabulatedRatingCurve / time' into the parameters" diff --git a/core/src/graph.jl b/core/src/graph.jl index 19215b250..edcd8f3cd 100644 --- a/core/src/graph.jl +++ b/core/src/graph.jl @@ -33,16 +33,15 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra node_id = NodeID(row.node_type, row.node_id) # Process allocation network ID if ismissing(row.subnetwork_id) - allocation_network_id = 0 + subnetwork_id = 0 else - allocation_network_id = row.subnetwork_id - if !haskey(node_ids, allocation_network_id) - node_ids[allocation_network_id] = Set{NodeID}() + subnetwork_id = row.subnetwork_id + if !haskey(node_ids, subnetwork_id) + node_ids[subnetwork_id] = Set{NodeID}() end - push!(node_ids[allocation_network_id], node_id) + push!(node_ids[subnetwork_id], node_id) end - graph[node_id] = - NodeMetadata(Symbol(snake_case(row.node_type)), allocation_network_id) + graph[node_id] = NodeMetadata(Symbol(snake_case(row.node_type)), subnetwork_id) end errors = false @@ -178,9 +177,16 @@ end """ Get the flow over the given edge (val is needed for get_tmp from ForwardDiff.jl). """ -function get_flow(graph::MetaGraph, id_src::NodeID, id_dst::NodeID, val)::Number - (; flow_dict, flow) = graph[] - return get_tmp(flow, val)[flow_dict[id_src, id_dst]] +function get_flow( + graph::MetaGraph, + id_src::NodeID, + id_dst::NodeID, + val; + prev::Bool = false, +)::Number + (; flow_dict, flow, flow_prev) = graph[] + flow_vector = prev ? flow_prev : flow + return get_tmp(flow_vector, val)[flow_dict[id_src, id_dst]] end """ diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 2566eabb4..c8de50fb3 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -43,7 +43,7 @@ const VectorInterpolation = """ Store information for a subnetwork used for allocation. -allocation_network_id: The ID of this allocation network +subnetwork_id: The ID of this allocation network capacity: The capacity per edge of the allocation network, as constrained by nodes that have a max_flow_rate problem: The JuMP.jl model for solving the allocation problem Δt_allocation: The time interval between consecutive allocation solves @@ -57,13 +57,15 @@ end """ Object for all information about allocation -allocation_network_ids: The unique sorted allocation network IDs +subnetwork_ids: The unique sorted allocation network IDs allocation models: The allocation models for the main network and subnetworks corresponding to - allocation_network_ids + subnetwork_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_demand: A record of demands and allocated flows for nodes that have these. +subnetwork_allocateds: The allocated flow of an edge from the main network to a subnetwork +mean_source_flows: Flows averaged over Δt_allocation over edges that are allocation sources +record_demand: A record of demands and allocated flows for nodes that have these record_flow: A record of all flows computed by allocation optimization, eventually saved to output file """ @@ -74,6 +76,7 @@ struct Allocation priorities::Vector{Int32} subnetwork_demands::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} subnetwork_allocateds::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} + mean_source_flows::Dict{Tuple{NodeID, NodeID}, Base.RefValue{Float64}} record_demand::@NamedTuple{ time::Vector{Float64}, subnetwork_id::Vector{Int32}, @@ -103,7 +106,7 @@ is_active(allocation::Allocation) = !isempty(allocation.allocation_models) """ Type for storing metadata of nodes in the graph type: type of the node -allocation_network_id: Allocation network ID (0 if not in subnetwork) +subnetwork_id: Allocation network ID (0 if not in subnetwork) """ struct NodeMetadata type::Symbol diff --git a/core/src/read.jl b/core/src/read.jl index 7d3ff0a69..060e3fa62 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -951,7 +951,7 @@ function Subgrid(db::DB, config::Config, basin::Basin)::Subgrid return Subgrid(basin_ids, interpolations, fill(NaN, length(basin_ids))) end -function Allocation(db::DB, config::Config)::Allocation +function Allocation(db::DB, config::Config, graph::MetaGraph)::Allocation record_demand = ( time = Float64[], subnetwork_id = Int32[], @@ -976,18 +976,27 @@ function Allocation(db::DB, config::Config)::Allocation optimization_type = String[], ) - allocation = Allocation( + mean_source_flows = Dict{Tuple{NodeID, NodeID}, Base.RefValue{Float64}}() + + # Find edges which serve as sources in allocation + for edge_metadata in values(graph.edge_data) + (; subnetwork_id_source, edge) = edge_metadata + if subnetwork_id_source != 0 + mean_source_flows[edge] = Ref(0.0) + end + end + + return Allocation( Int32[], AllocationModel[], Vector{Tuple{NodeID, NodeID}}[], get_all_priorities(db, config), - Dict{Tuple{NodeID, NodeID}, Float64}(), - Dict{Tuple{NodeID, NodeID}, Float64}(), + Dict{Tuple{NodeID, NodeID}, Vector{Float64}}(), + Dict{Tuple{NodeID, NodeID}, Vector{Float64}}(), + mean_source_flows, record_demand, record_flow, ) - - return allocation end """ @@ -1007,7 +1016,7 @@ function Parameters(db::DB, config::Config)::Parameters n_states = length(get_ids(db, "Basin")) + length(get_ids(db, "PidControl")) chunk_sizes = get_chunk_sizes(config, n_states) graph = create_graph(db, config, chunk_sizes) - allocation = Allocation(db, config) + allocation = Allocation(db, config, graph) if !valid_edges(graph) error("Invalid edge(s) found.") From 86a87ffcfdc40c93eebc78f25519e80060694903 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 24 Apr 2024 11:22:51 +0200 Subject: [PATCH 2/4] Fix tests --- core/test/allocation_test.jl | 34 ++++++++++++++++++++++------------ 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 023f9dee4..1ad3c8816 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -11,10 +11,10 @@ db = SQLite.DB(db_path) p = Ribasim.Parameters(db, cfg) - graph = p.graph + (; graph, allocation) = p close(db) - Ribasim.set_flow!(graph, NodeID(:FlowBoundary, 1), NodeID(:Basin, 2), 4.5) # Source flow + allocation.mean_source_flows[(NodeID(:FlowBoundary, 1), NodeID(:Basin, 2))][] = 4.5 allocation_model = p.allocation.allocation_models[1] u = ComponentVector(; storage = zeros(length(p.basin.node_id))) Ribasim.allocate!(p, allocation_model, 0.0, u, OptimizationType.allocate) @@ -210,8 +210,13 @@ end close(db) (; allocation, user_demand, graph, basin) = p - (; allocation_models, subnetwork_demands, subnetwork_allocateds, record_flow) = - allocation + (; + allocation_models, + subnetwork_demands, + subnetwork_allocateds, + record_flow, + mean_source_flows, + ) = allocation t = 0.0 # Collecting demands @@ -241,7 +246,8 @@ end @test F_abs_user_demand[NodeID(:Pump, 38)] ∈ objective_variables # Running full allocation algorithm - Ribasim.set_flow!(graph, NodeID(:FlowBoundary, 1), NodeID(:Basin, 2), 4.5) + (; Δt_allocation) = allocation_models[1] + mean_source_flows[(NodeID(:FlowBoundary, 1), NodeID(:Basin, 2))][] = 4.5 * Δt_allocation u = ComponentVector(; storage = zeros(length(p.basin.node_id))) Ribasim.update_allocation!((; p, t, u)) @@ -283,12 +289,13 @@ end close(db) (; allocation, user_demand, graph, basin) = p - (; allocation_models, subnetwork_demands, subnetwork_allocateds) = allocation + (; allocation_models, subnetwork_demands, subnetwork_allocateds, mean_source_flows) = + allocation t = 0.0 - # Set flows of sources in subnetworks - Ribasim.set_flow!(graph, NodeID(:FlowBoundary, 58), NodeID(:Basin, 16), 1.0) - Ribasim.set_flow!(graph, NodeID(:FlowBoundary, 59), NodeID(:Basin, 44), 1e-3) + # Set flows of sources in + mean_source_flows[(NodeID(:FlowBoundary, 58), NodeID(:Basin, 16))][] = 1.0 + mean_source_flows[(NodeID(:FlowBoundary, 59), NodeID(:Basin, 44))][] = 1e-3 # Collecting demands u = ComponentVector(; storage = zeros(length(basin.node_id))) @@ -372,7 +379,7 @@ end @test storage[stage_6] ≈ u_stage_6.(t[stage_6]) rtol = 1e-4 end -@testitem "flow_demand" begin +@testitem "Flow demand" begin using JuMP using Ribasim: NodeID, OptimizationType @@ -422,6 +429,9 @@ end t = 0.0 (; u) = model.integrator optimization_type = OptimizationType.internal_sources + for (edge, value) in allocation.mean_source_flows + value[] = Ribasim.get_flow(graph, edge..., 0) + end Ribasim.set_initial_values!(allocation_model, p, u, t) # Priority 1 @@ -484,8 +494,8 @@ end ) # Get demand from buffers d = user_demand.demand_itp[3][4](t) - @assert JuMP.value(F[(NodeID(:UserDemand, 4), NodeID(:Basin, 7))]) + - JuMP.value(F[(NodeID(:UserDemand, 6), NodeID(:Basin, 7))]) == d + @test JuMP.value(F[(NodeID(:UserDemand, 4), NodeID(:Basin, 7))]) + + JuMP.value(F[(NodeID(:UserDemand, 6), NodeID(:Basin, 7))]) == d end @testitem "flow_demand_with_max_flow_rate" begin From 675a0e41053950a33614cfdb26b70ce6689c13ee Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 24 Apr 2024 15:17:52 +0200 Subject: [PATCH 3/4] Use average vertical fluxes on basin for level demand --- core/src/allocation_init.jl | 8 +++++--- core/src/allocation_optim.jl | 20 ++++++++++++-------- core/src/callback.jl | 26 ++++++++++++++++++-------- core/src/parameter.jl | 4 ++-- core/src/read.jl | 14 +++++++++++--- core/src/util.jl | 7 ++++--- core/test/allocation_test.jl | 14 +++++++------- 7 files changed, 59 insertions(+), 34 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 9bf3ddb30..eea232e5c 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -153,7 +153,8 @@ function add_variables_basin!( (; graph) = p node_ids_basin = [ node_id for - node_id in graph[].node_ids[subnetwork_id] if graph[node_id].type == :basin + node_id in graph[].node_ids[subnetwork_id] if graph[node_id].type == :basin && + has_external_demand(graph, node_id, :level_demand)[1] ] problem[:F_basin_in] = JuMP.@variable(problem, F_basin_in[node_id = node_ids_basin,] >= 0.0) @@ -211,7 +212,7 @@ function add_variables_absolute_value!( type = node_id.type if type == NodeType.UserDemand push!(node_ids_user_demand, node_id) - elseif type == NodeType.Basin + elseif has_external_demand(graph, node_id, :level_demand)[1] push!(node_ids_level_demand, node_id) elseif has_external_demand(graph, node_id, :flow_demand)[1] push!(node_ids_flow_demand, node_id) @@ -393,7 +394,8 @@ function add_constraints_conservation_node!( end # If the node is a Basin, add basin in- and outflow - if node_id.type == NodeType.Basin + if node_id.type == NodeType.Basin && + has_external_demand(graph, node_id, :level_demand)[1] push!(inflows_node, F_basin_out[node_id]) push!(outflows_node, F_basin_in[node_id]) end diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 1653c989e..2103f2b9f 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -255,7 +255,7 @@ function set_initial_capacities_source!( )::Nothing (; problem) = allocation_model (; graph, allocation) = p - (; mean_source_flows) = allocation + (; mean_flows) = allocation (; subnetwork_id) = allocation_model source_constraints = problem[:source] main_network_source_edges = get_main_network_connections(p, subnetwork_id) @@ -266,7 +266,7 @@ function set_initial_capacities_source!( # If it is a source edge for this allocation problem if edge ∉ main_network_source_edges # Reset the source to the averaged flow over the last allocation period - source_capacity = mean_source_flows[edge][] + source_capacity = mean_flows[edge][] JuMP.set_normalized_rhs( source_constraints[edge], # It is assumed that the allocation procedure does not have to be differentiated. @@ -358,14 +358,14 @@ function get_basin_data( u::ComponentVector, node_id::NodeID, ) - (; graph, basin, level_demand) = p + (; graph, basin, level_demand, allocation) = p (; vertical_flux) = basin (; Δt_allocation) = allocation_model + (; mean_flows) = allocation @assert node_id.type == NodeType.Basin vertical_flux = get_tmp(vertical_flux, 0) _, basin_idx = id_index(basin.node_id, node_id) - # NOTE: Instantaneous - influx = get_influx(basin, node_id) + influx = mean_flows[(node_id, node_id)][] _, 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) @@ -512,12 +512,15 @@ function set_initial_demands_level!( p::Parameters, t::Float64, )::Nothing - (; subnetwork_id) = allocation_model + (; subnetwork_id, problem) = allocation_model (; graph, basin) = p (; node_id, demand) = basin - for (i, id) in enumerate(node_id) + node_ids_level_demand = only(problem[:basin_outflow].axes) + + for id in node_ids_level_demand if graph[id].subnetwork_id == subnetwork_id + _, i = id_index(node_id, id) demand[i] = get_basin_demand(allocation_model, u, p, t, id) end end @@ -607,8 +610,9 @@ function adjust_demands_level!(allocation_model::AllocationModel, p::Parameters) F_basin_in = problem[:F_basin_in] # Reduce the demand by what was allocated - for (i, id) in enumerate(node_id) + for id in only(F_basin_in.axes) if graph[id].subnetwork_id == subnetwork_id + _, i = id_index(basin.node_id, id) demand[i] -= JuMP.value(F_basin_in[id]) end end diff --git a/core/src/callback.jl b/core/src/callback.jl index b22b8901e..048ba7530 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -114,11 +114,21 @@ function integrate_flows!(u, t, integrator)::Nothing end # Allocation source flows - for (edge, value) in allocation.mean_source_flows - value[] += - 0.5 * - (get_flow(graph, edge..., 0) + get_flow(graph, edge..., 0; prev = true)) * - dt + for (edge, value) in allocation.mean_flows + if edge[1] == edge[2] + # Vertical fluxes + _, basin_idx = id_index(basin.node_id, edge[1]) + value[] += + 0.5 * + (get_influx(basin, basin_idx) + get_influx(basin, basin_idx; prev = true)) * + dt + else + # Horizontal flows + value[] += + 0.5 * + (get_flow(graph, edge..., 0) + get_flow(graph, edge..., 0; prev = true)) * + dt + end end copyto!(flow_prev, flow) @@ -455,12 +465,12 @@ end function update_allocation!(integrator)::Nothing (; p, t, u) = integrator (; allocation) = p - (; allocation_models, mean_source_flows) = allocation + (; allocation_models, mean_flows) = allocation (; Δt_allocation) = allocation_models[1] # Divide by the allocation Δt to obtain the mean flows # from the integrated flows - for value in values(mean_source_flows) + for value in values(mean_flows) value[] /= Δt_allocation end @@ -480,7 +490,7 @@ function update_allocation!(integrator)::Nothing end # Reset the mean source flows - for value in values(mean_source_flows) + for value in values(mean_flows) value[] = 0.0 end end diff --git a/core/src/parameter.jl b/core/src/parameter.jl index c8de50fb3..d16b793a8 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -64,7 +64,7 @@ main_network_connections: (from_id, to_id) from the main network to the subnetwo priorities: All used priority values. subnetwork_demands: The demand of an edge from the main network to a subnetwork subnetwork_allocateds: The allocated flow of an edge from the main network to a subnetwork -mean_source_flows: Flows averaged over Δt_allocation over edges that are allocation sources +mean_flows: Flows averaged over Δt_allocation over edges that are allocation sources record_demand: A record of demands and allocated flows for nodes that have these record_flow: A record of all flows computed by allocation optimization, eventually saved to output file @@ -76,7 +76,7 @@ struct Allocation priorities::Vector{Int32} subnetwork_demands::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} subnetwork_allocateds::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} - mean_source_flows::Dict{Tuple{NodeID, NodeID}, Base.RefValue{Float64}} + mean_flows::Dict{Tuple{NodeID, NodeID}, Base.RefValue{Float64}} record_demand::@NamedTuple{ time::Vector{Float64}, subnetwork_id::Vector{Int32}, diff --git a/core/src/read.jl b/core/src/read.jl index 060e3fa62..f9470f7c0 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -976,13 +976,21 @@ function Allocation(db::DB, config::Config, graph::MetaGraph)::Allocation optimization_type = String[], ) - mean_source_flows = Dict{Tuple{NodeID, NodeID}, Base.RefValue{Float64}}() + mean_flows = Dict{Tuple{NodeID, NodeID}, Base.RefValue{Float64}}() # Find edges which serve as sources in allocation for edge_metadata in values(graph.edge_data) (; subnetwork_id_source, edge) = edge_metadata if subnetwork_id_source != 0 - mean_source_flows[edge] = Ref(0.0) + mean_flows[edge] = Ref(0.0) + end + end + + # Find basins with a level demand + for node_id in values(graph.vertex_labels) + if node_id.type == NodeType.Basin && + has_external_demand(graph, node_id, :level_demand)[1] + mean_flows[(node_id, node_id)] = Ref(0.0) end end @@ -993,7 +1001,7 @@ function Allocation(db::DB, config::Config, graph::MetaGraph)::Allocation get_all_priorities(db, config), Dict{Tuple{NodeID, NodeID}, Vector{Float64}}(), Dict{Tuple{NodeID, NodeID}, Vector{Float64}}(), - mean_source_flows, + mean_flows, record_demand, record_flow, ) diff --git a/core/src/util.jl b/core/src/util.jl index a52bdb277..b1f816d9c 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -711,10 +711,11 @@ function get_influx(basin::Basin, node_id::NodeID)::Float64 return get_influx(basin, basin_idx) end -function get_influx(basin::Basin, basin_idx::Int)::Float64 - (; vertical_flux) = basin +function get_influx(basin::Basin, basin_idx::Int; prev::Bool = false)::Float64 + (; vertical_flux, vertical_flux_prev) = basin vertical_flux = get_tmp(vertical_flux, 0) - (; precipitation, evaporation, drainage, infiltration) = vertical_flux + flux_vector = prev ? vertical_flux_prev : vertical_flux + (; precipitation, evaporation, drainage, infiltration) = flux_vector return precipitation[basin_idx] - evaporation[basin_idx] + drainage[basin_idx] - infiltration[basin_idx] end diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 0a5069cc8..d49f0512d 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -14,7 +14,7 @@ (; graph, allocation) = p close(db) - allocation.mean_source_flows[(NodeID(:FlowBoundary, 1), NodeID(:Basin, 2))][] = 4.5 + allocation.mean_flows[(NodeID(:FlowBoundary, 1), NodeID(:Basin, 2))][] = 4.5 allocation_model = p.allocation.allocation_models[1] u = ComponentVector(; storage = zeros(length(p.basin.node_id))) Ribasim.allocate!(p, allocation_model, 0.0, u, OptimizationType.allocate) @@ -215,7 +215,7 @@ end subnetwork_demands, subnetwork_allocateds, record_flow, - mean_source_flows, + mean_flows, ) = allocation t = 0.0 @@ -247,7 +247,7 @@ end # Running full allocation algorithm (; Δt_allocation) = allocation_models[1] - mean_source_flows[(NodeID(:FlowBoundary, 1), NodeID(:Basin, 2))][] = 4.5 * Δt_allocation + mean_flows[(NodeID(:FlowBoundary, 1), NodeID(:Basin, 2))][] = 4.5 * Δt_allocation u = ComponentVector(; storage = zeros(length(p.basin.node_id))) Ribasim.update_allocation!((; p, t, u)) @@ -289,13 +289,13 @@ end close(db) (; allocation, user_demand, graph, basin) = p - (; allocation_models, subnetwork_demands, subnetwork_allocateds, mean_source_flows) = + (; allocation_models, subnetwork_demands, subnetwork_allocateds, mean_flows) = allocation t = 0.0 # Set flows of sources in - mean_source_flows[(NodeID(:FlowBoundary, 58), NodeID(:Basin, 16))][] = 1.0 - mean_source_flows[(NodeID(:FlowBoundary, 59), NodeID(:Basin, 44))][] = 1e-3 + mean_flows[(NodeID(:FlowBoundary, 58), NodeID(:Basin, 16))][] = 1.0 + mean_flows[(NodeID(:FlowBoundary, 59), NodeID(:Basin, 44))][] = 1e-3 # Collecting demands u = ComponentVector(; storage = zeros(length(basin.node_id))) @@ -429,7 +429,7 @@ end t = 0.0 (; u) = model.integrator optimization_type = OptimizationType.internal_sources - for (edge, value) in allocation.mean_source_flows + for (edge, value) in allocation.mean_flows value[] = Ribasim.get_flow(graph, edge..., 0) end Ribasim.set_initial_values!(allocation_model, p, u, t) From 2dbf9cfaa56d4b941561b3878f1a00be440933cd Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 24 Apr 2024 16:26:46 +0200 Subject: [PATCH 4/4] Update documentation --- core/src/allocation_init.jl | 5 ++--- core/src/allocation_optim.jl | 3 +-- core/src/read.jl | 3 +-- docs/core/allocation.qmd | 14 +++++++++----- docs/core/equations.qmd | 2 +- 5 files changed, 14 insertions(+), 13 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index eea232e5c..6ac9f2be0 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -393,9 +393,8 @@ function add_constraints_conservation_node!( end end - # If the node is a Basin, add basin in- and outflow - if node_id.type == NodeType.Basin && - has_external_demand(graph, node_id, :level_demand)[1] + # If the node is a Basin with a level demand, add basin in- and outflow + if has_external_demand(graph, node_id, :level_demand)[1] push!(inflows_node, F_basin_out[node_id]) push!(outflows_node, F_basin_in[node_id]) end diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 2103f2b9f..01370c228 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -768,8 +768,7 @@ function save_demands_and_allocations!( #NOTE: instantaneous realized = get_flow(graph, inflow_id(graph, node_id), node_id, 0) - elseif node_id.type == NodeType.Basin && - has_external_demand(graph, node_id, :level_demand)[1] + elseif has_external_demand(graph, node_id, :level_demand)[1] basin_priority_idx = get_external_priority_idx(p, node_id) if priority_idx == 1 || basin_priority_idx == priority_idx diff --git a/core/src/read.jl b/core/src/read.jl index f9470f7c0..9b84fed30 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -988,8 +988,7 @@ function Allocation(db::DB, config::Config, graph::MetaGraph)::Allocation # Find basins with a level demand for node_id in values(graph.vertex_labels) - if node_id.type == NodeType.Basin && - has_external_demand(graph, node_id, :level_demand)[1] + if has_external_demand(graph, node_id, :level_demand)[1] mean_flows[(node_id, node_id)] = Ref(0.0) end end diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index 9a5947482..1399513b3 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -42,7 +42,11 @@ Sources are indicated by a set of edges in the subnetwork $$ E_S^\text{source} \subset E. $$ -That is, if $(i,j) \in E_S^\text{source}$, then $Q_{ij}$ (see the [formal model description](equations.qmd#formal-model-description)) is treated as a source flow in the allocation problem. These edges are either coming from a boundary/source node (e.g. a level or flow boundary) or connect the main network to a subnetwork. +That is, if $(i,j) \in E_S^\text{source}$, then the average over the allocation interval $\Delta t_{\text{alloc}}$ of the of the flow over this edge +$$ + \frac{1}{\Delta t_{\text{alloc}}}\int_{t - \Delta t_{\text{alloc}}}^tQ_{ij}(t') dt' +$$ + is treated as a source flow in the allocation problem. These edges are either coming from a boundary/source node (e.g. a level or flow boundary) or connect the main network to a subnetwork. For the definition of $Q_{ij}$ see [the formal model description](equations.qmd#formal-model-description). ### User demands @@ -71,12 +75,12 @@ for all $i \in FD_S$. Here $d^{p_{\text{df}}}$ is given by the original flow dem ### Vertical fluxes and local storage -Apart from the source flows denoted by edges, there are other sources of water in the subnetwork, associated with the basins in the subnetwork $B_S = B \cap S$. Firstly there is the sum of the vertical fluxes (precipitation, evaporation, infiltration and drainage) for each basin +Apart from the source flows denoted by edges, there are other sources of water in the subnetwork, associated with the basins in the subnetwork $B_S = B \cap S$. Firstly there is the average over the allocation interval $\Delta t_{\text{alloc}}$ of the vertical fluxes (precipitation, evaporation, infiltration and drainage) for each basin: $$ - \phi_i(t), \quad \forall i \in B_S. + \phi_i(t) = \frac{1}{\Delta t_{\text{alloc}}}\int_{t - \Delta t_{\text{alloc}}}^t \left[Q_{P,i}(t') - Q_{E,i}(t') + Q_{\text{drn},i}(t') - Q_{\text{inf},i}(t') \right] dt', \quad \forall i \in B_S. $$ -We consider fluxes into the basin to be positive and out of the basin to be negative. +We consider fluxes into the basin to be positive and out of the basin to be negative. For more information see [the natural water balance terms](equations.qmd#natural-water-balance-terms). 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 $$ @@ -88,7 +92,7 @@ $$ 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$. +for all $i \in B_S$. Note that the basin demand has only a single priority, so for other priorities this demand is $0$. ### Constraining factors diff --git a/docs/core/equations.qmd b/docs/core/equations.qmd index 1f395cd58..a2719a638 100644 --- a/docs/core/equations.qmd +++ b/docs/core/equations.qmd @@ -137,7 +137,7 @@ ax.hlines([0,1],x_min,x_max, color = "k", ls = ":", zorder=-1) ax.vlines([0,p], *y_lim, color = "k", ls = ":") ax.set_xlim(x_min,x_max) ax.set_xlabel("$x$", fontsize=fontsize) -ax.set_ylabel("$\phi(x;p)$", fontsize=fontsize) +ax.set_ylabel(r"$\phi(x;p)$", fontsize=fontsize) ax.set_ylim(y_lim) fig.tight_layout()