From 26d90a689afed217a155f00ce02d737148417683 Mon Sep 17 00:00:00 2001 From: Bart de Koning <74617371+SouthEndMusic@users.noreply.github.com> Date: Thu, 19 Oct 2023 16:36:21 +0200 Subject: [PATCH] Allocation beyond max flow (#663) Fixes https://github.com/Deltares/Ribasim/issues/632, fixes https://github.com/Deltares/Ribasim/issues/14, fixes https://github.com/Deltares/Ribasim/issues/349, using `JuMP.jl` and `HiGHS`. --------- Co-authored-by: Martijn Visser --- core/Manifest.toml | 84 +- core/Project.toml | 4 + core/src/Ribasim.jl | 17 +- core/src/allocation.jl | 833 ++++++++++++++++++ core/src/create.jl | 4 + core/src/solve.jl | 27 + core/src/utils.jl | 45 + core/test/allocation.jl | 34 + core/test/runtests.jl | 1 + core/test/time.jl | 1 - docs/_quarto.yml | 1 + docs/core/allocation.qmd | 223 +++++ docs/core/index.qmd | 2 +- python/ribasim/ribasim/node_types/user.py | 4 +- .../ribasim_testmodels/__init__.py | 7 +- .../ribasim_testmodels/allocation.py | 264 +++++- .../ribasim_testmodels/dutch_waterways.py | 13 - 17 files changed, 1537 insertions(+), 27 deletions(-) create mode 100644 core/src/allocation.jl create mode 100644 core/test/allocation.jl create mode 100644 docs/core/allocation.qmd diff --git a/core/Manifest.toml b/core/Manifest.toml index d04de9383..fd14b4ebc 100644 --- a/core/Manifest.toml +++ b/core/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.9.3" manifest_format = "2.0" -project_hash = "40e443b575cd5dffce730be3f39503b3b068d36a" +project_hash = "5ce3afdbcba98d75c8a7887863cdc1ee1cabf1e7" [[deps.ADTypes]] git-tree-sha1 = "5d2e21d7b0d8c22f67483ef95ebdc39c0e6b6003" @@ -85,6 +85,12 @@ git-tree-sha1 = "0e2855d28cc3983a9edf4bee18478b020178a43f" uuid = "59605e27-edc0-445a-b93d-c09a3a50b330" version = "0.1.0" +[[deps.BenchmarkTools]] +deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"] +git-tree-sha1 = "d9a9701b899b30332bbcb3e1679c41cce81fb0e8" +uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +version = "1.3.2" + [[deps.BitIntegers]] deps = ["Random"] git-tree-sha1 = "a55462dfddabc34bc97d3a7403a2ca2802179ae6" @@ -97,6 +103,12 @@ git-tree-sha1 = "0c5f81f47bbbcf4aea7b2959135713459170798b" uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b" version = "0.1.5" +[[deps.Bzip2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2" +uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" +version = "1.0.8+0" + [[deps.CEnum]] git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90" uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" @@ -124,12 +136,24 @@ git-tree-sha1 = "70232f82ffaab9dc52585e0dd043b5e0c6b714f1" uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9" version = "0.1.12" +[[deps.CodecBzip2]] +deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"] +git-tree-sha1 = "c0ae2a86b162fb5d7acc65269b469ff5b8a73594" +uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd" +version = "0.8.1" + [[deps.CodecLz4]] deps = ["Lz4_jll", "TranscodingStreams"] git-tree-sha1 = "8bf4f9e2ee52b5e217451a7cd9171fcd4e16ae23" uuid = "5ba52731-8f18-5e0d-9241-30f10d1ec561" version = "0.4.1" +[[deps.CodecZlib]] +deps = ["TranscodingStreams", "Zlib_jll"] +git-tree-sha1 = "cd67fc487743b2f0fd4380d4cbd3a24660d0eec8" +uuid = "944b1d66-785c-5afd-91f1-9de20f533193" +version = "0.7.3" + [[deps.CodecZstd]] deps = ["CEnum", "TranscodingStreams", "Zstd_jll"] git-tree-sha1 = "849470b337d0fa8449c21061de922386f32949d9" @@ -475,6 +499,18 @@ git-tree-sha1 = "899050ace26649433ef1af25bc17a815b3db52b7" uuid = "86223c79-3864-5bf0-83f7-82e725a168b6" version = "1.9.0" +[[deps.HiGHS]] +deps = ["HiGHS_jll", "MathOptInterface", "PrecompileTools", "SparseArrays"] +git-tree-sha1 = "9d75ef949c17a2a150b91b8365a6e5bc43a2a0d3" +uuid = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" +version = "1.7.3" + +[[deps.HiGHS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] +git-tree-sha1 = "10bf0ecdf70f643bfc1948a6af0a98be3950a3fc" +uuid = "8fd58aa0-07eb-5a78-9b36-339c94fd15ea" +version = "1.6.0+0" + [[deps.HostCPUFeatures]] deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"] git-tree-sha1 = "eb8fed28f4994600e29beef49744639d985a04b2" @@ -533,6 +569,24 @@ git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca" uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" version = "1.5.0" +[[deps.JSON]] +deps = ["Dates", "Mmap", "Parsers", "Unicode"] +git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a" +uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +version = "0.21.4" + +[[deps.JuMP]] +deps = ["LinearAlgebra", "MacroTools", "MathOptInterface", "MutableArithmetics", "OrderedCollections", "Printf", "SnoopPrecompile", "SparseArrays"] +git-tree-sha1 = "3700a700bc80856fe673b355123ae4574f2d5dfe" +uuid = "4076af6c-e467-56ae-b986-b466b2749572" +version = "1.15.1" + + [deps.JuMP.extensions] + JuMPDimensionalDataExt = "DimensionalData" + + [deps.JuMP.weakdeps] + DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0" + [[deps.KLU]] deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse_jll"] git-tree-sha1 = "884c2968c2e8e7e6bf5956af88cb46aa745c854b" @@ -700,6 +754,12 @@ version = "0.1.8" deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" +[[deps.MathOptInterface]] +deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test", "Unicode"] +git-tree-sha1 = "5c9f1e635e8d491297e596b56fec1c95eafb95a3" +uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +version = "1.20.1" + [[deps.MbedTLS_jll]] deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" @@ -723,6 +783,12 @@ git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab" uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" version = "0.2.4" +[[deps.MutableArithmetics]] +deps = ["LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "6985021d02ab8c509c841bb8b2becd3145a7b490" +uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +version = "1.3.3" + [[deps.NLSolversBase]] deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"] git-tree-sha1 = "a0b464d183da839699f4c79e7606d9d186ec172c" @@ -747,9 +813,17 @@ version = "1.2.0" [[deps.NonlinearSolve]] deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "EnumX", "FiniteDiff", "ForwardDiff", "LineSearches", "LinearAlgebra", "LinearSolve", "PrecompileTools", "RecursiveArrayTools", "Reexport", "SciMLBase", "SimpleNonlinearSolve", "SparseArrays", "SparseDiffTools", "StaticArraysCore", "UnPack"] -git-tree-sha1 = "a5f1f836da05d513c4143576af8f5d8e51b759f5" +git-tree-sha1 = "ee92770e0832314ccd424d83a0ab4c75fc6dc91f" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -version = "2.2.1" +version = "2.3.0" + + [deps.NonlinearSolve.extensions] + NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt" + NonlinearSolveLeastSquaresOptimExt = "LeastSquaresOptim" + + [deps.NonlinearSolve.weakdeps] + FastLevenbergMarquardt = "7a0df574-e128-4d35-8cbd-3d84502bf7ce" + LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891" [[deps.OffsetArrays]] deps = ["Adapt"] @@ -853,6 +927,10 @@ version = "1.4.1" deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" +[[deps.Profile]] +deps = ["Printf"] +uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79" + [[deps.ProgressLogging]] deps = ["Logging", "SHA", "UUIDs"] git-tree-sha1 = "80d919dee55b9c50e8d9e2da5eeafff3fe58b539" diff --git a/core/Project.toml b/core/Project.toml index 31cf8981c..684f36f6b 100644 --- a/core/Project.toml +++ b/core/Project.toml @@ -19,7 +19,9 @@ DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" Legolas = "741b9549-f6ed-4911-9fbf-4a1c0c97f0cd" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36" @@ -50,7 +52,9 @@ DiffEqCallbacks = "2.29.1" FiniteDiff = "2.21" ForwardDiff = "0.10" Graphs = "1.6" +HiGHS = "1.7" IterTools = "1.4" +JuMP = "1.15" Legolas = "0.5" LoggingExtras = "1" OrdinaryDiffEq = "6.7" diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index 1acfde3f4..6041bb6f3 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -16,6 +16,8 @@ module Ribasim import IterTools import BasicModelInterface as BMI +import HiGHS +import JuMP.Model as JuMPModel import TranscodingStreams using Arrow: Arrow, Table @@ -29,7 +31,19 @@ using DBInterface: execute, prepare using Dictionaries: Indices, Dictionary, gettoken, dictionary using ForwardDiff: pickchunksize using DiffEqCallbacks -using Graphs: DiGraph, add_edge!, adjacency_matrix, inneighbors, outneighbors +using Graphs: + add_edge!, + adjacency_matrix, + all_neighbors, + DiGraph, + Edge, + edges, + inneighbors, + nv, + outneighbors, + rem_edge! + +using JuMP: @variable, @constraint, @objective, set_normalized_rhs, optimize!, value using Legolas: Legolas, @schema, @version, validate, SchemaVersion, declared using Logging: current_logger, min_enabled_level, with_logger using LoggingExtras: EarlyFilteredLogger, LevelOverrideLogger @@ -48,6 +62,7 @@ TimerOutputs.complement!() include("validation.jl") include("solve.jl") +include("allocation.jl") include("config.jl") using .config include("utils.jl") diff --git a/core/src/allocation.jl b/core/src/allocation.jl new file mode 100644 index 000000000..bef4f28d8 --- /dev/null +++ b/core/src/allocation.jl @@ -0,0 +1,833 @@ +""" +Get: +- The mapping from subnetwork node IDs to allocation graph node IDs +- The mapping from allocation graph source node IDs to subnetwork source edge IDs +""" +function get_node_id_mapping( + p::Parameters, + subnetwork_node_ids::Vector{Int}, + source_edge_ids::Vector{Int}, +) + (; lookup, connectivity) = p + (; graph_flow, edge_ids_flow_inv) = connectivity + + # Mapping node_id => (allocgraph_node_id, type) where such a correspondence exists; + # allocgraph_node_type in [:user, :junction, :basin, :source] + node_id_mapping = Dict{Int, Tuple{Int, Symbol}}() + + # Determine the number of nodes in the allocgraph + n_allocgraph_nodes = 0 + for subnetwork_node_id in subnetwork_node_ids + add_allocgraph_node = false + node_type = lookup[subnetwork_node_id] + + if node_type in [:user, :basin] + add_allocgraph_node = true + + elseif length(all_neighbors(graph_flow, subnetwork_node_id)) > 2 + # Each junction (that is, a node with more than 2 neighbors) + # in the subnetwork gets an allocgraph node + add_allocgraph_node = true + node_type = :junction + end + + if add_allocgraph_node + n_allocgraph_nodes += 1 + node_id_mapping[subnetwork_node_id] = (n_allocgraph_nodes, node_type) + end + end + + # Add nodes in the allocation graph for nodes connected in the problem to the source edges + # One of these nodes can be outside the subnetwork, as long as the edge + # connects to the subnetwork + # Source edge mapping: allocation graph source node ID => subnetwork source edge ID + source_edge_mapping = Dict{Int, Int}() + for source_edge_id in source_edge_ids + subnetwork_node_id_1, subnetwork_node_id_2 = edge_ids_flow_inv[source_edge_id] + if subnetwork_node_id_1 ∉ keys(node_id_mapping) + n_allocgraph_nodes += 1 + node_id_mapping[subnetwork_node_id_1] = (n_allocgraph_nodes, :source) + source_edge_mapping[n_allocgraph_nodes] = source_edge_id + else + node_id_mapping[subnetwork_node_id_1][2] = :source + source_edge_mapping[n_allocgraph_nodes] = source_edge_id + end + if subnetwork_node_id_2 ∉ keys(node_id_mapping) + n_allocgraph_nodes += 1 + node_id_mapping[subnetwork_node_id_2] = (n_allocgraph_nodes, :junction) + end + end + return node_id_mapping, source_edge_mapping +end + +""" +This loop finds allocgraph edges in several ways: +- Between allocgraph nodes whose equivalent in the subnetwork are directly connected +- Between allocgraph nodes whose equivalent in the subnetwork are connected + with one or more non-junction nodes in between + +Here edges are added to the allocation graph that are given by a single edge in +the subnetwork. +""" +function find_allocation_graph_edges!( + graph_allocation::DiGraph{Int}, + node_id_mapping::Dict{Int, Tuple{Int, Symbol}}, + p::Parameters, + subnetwork_node_ids::Vector{Int}, +)::Tuple{Vector{Vector{Int}}, SparseMatrixCSC{Float64, Int}} + (; connectivity, user) = p + (; graph_flow) = connectivity + + allocgraph_edges_composite = Vector{Int}[] + n_allocgraph_nodes = nv(graph_allocation) + capacity = spzeros(n_allocgraph_nodes, n_allocgraph_nodes) + + for subnetwork_node_id in subnetwork_node_ids + subnetwork_inneighbor_ids = inneighbors(graph_flow, subnetwork_node_id) + subnetwork_outneighbor_ids = outneighbors(graph_flow, subnetwork_node_id) + subnetwork_neighbor_ids = all_neighbors(graph_flow, subnetwork_node_id) + + if subnetwork_node_id in keys(node_id_mapping) + if subnetwork_node_id ∉ user.node_id + # Direct connections in the subnetwork between nodes that + # have an equivalent allocgraph graph node + for subnetwork_inneighbor_id in subnetwork_inneighbor_ids + if subnetwork_inneighbor_id in keys(node_id_mapping) + allocgraph_node_id_1 = node_id_mapping[subnetwork_node_id][1] + allocgraph_node_id_2 = node_id_mapping[subnetwork_inneighbor_id][1] + add_edge!( + graph_allocation, + allocgraph_node_id_2, + allocgraph_node_id_1, + ) + # These direct connections cannot have capacity constraints + capacity[allocgraph_node_id_2, allocgraph_node_id_1] = Inf + end + end + for subnetwork_outneighbor_id in subnetwork_outneighbor_ids + if subnetwork_outneighbor_id in keys(node_id_mapping) + allocgraph_node_id_1 = node_id_mapping[subnetwork_node_id][1] + allocgraph_node_id_2 = node_id_mapping[subnetwork_outneighbor_id][1] + add_edge!( + graph_allocation, + allocgraph_node_id_1, + allocgraph_node_id_2, + ) + # if subnetwork_outneighbor_id in user.node_id: Capacity depends on user demand at a given priority + # else: These direct connections cannot have capacity constraints + capacity[allocgraph_node_id_1, allocgraph_node_id_2] = Inf + end + end + end + else + # Try to find an existing allocgraph composite edge to add the current subnetwork_node_id to + found_edge = false + for allocgraph_edge_composite in allocgraph_edges_composite + if allocgraph_edge_composite[1] in subnetwork_neighbor_ids + pushfirst!(allocgraph_edge_composite, subnetwork_node_id) + found_edge = true + break + elseif allocgraph_edge_composite[end] in subnetwork_neighbor_ids + push!(allocgraph_edge_composite, subnetwork_node_id) + found_edge = true + break + end + end + + # Start a new allocgraph composite edge if no existing edge to append to was found + if !found_edge + push!(allocgraph_edges_composite, [subnetwork_node_id]) + end + end + end + return allocgraph_edges_composite, capacity +end + +""" +For the composite allocgraph edges: +- Find out whether they are connected to allocgraph nodes on both ends +- Compute their capacity +- Find out their allowed flow direction(s) +""" +function process_allocation_graph_edges!( + graph_allocation::DiGraph{Int}, + capacity::SparseMatrixCSC{Float64, Int}, + allocgraph_edges_composite::Vector{Vector{Int}}, + node_id_mapping::Dict{Int, Tuple{Int, Symbol}}, + p::Parameters, +)::SparseMatrixCSC{Float64, Int} + (; connectivity, lookup) = p + (; graph_flow) = connectivity + n_allocgraph_nodes = nv(graph_allocation) + + for allocgraph_edge_composite in allocgraph_edges_composite + # Find allocgraph node connected to this edge on the first end + allocgraph_node_id_1 = nothing + subnetwork_neighbors_side_1 = + all_neighbors(graph_flow, allocgraph_edge_composite[1]) + for subnetwork_neighbor_node_id in subnetwork_neighbors_side_1 + if subnetwork_neighbor_node_id in keys(node_id_mapping) + allocgraph_node_id_1 = node_id_mapping[subnetwork_neighbor_node_id][1] + pushfirst!(allocgraph_edge_composite, subnetwork_neighbor_node_id) + break + end + end + + # No connection to a max flow node found on this side, so edge is discarded + if isnothing(allocgraph_node_id_1) + continue + end + + # Find allocgraph node connected to this edge on the second end + allocgraph_node_id_2 = nothing + subnetwork_neighbors_side_2 = + all_neighbors(graph_flow, allocgraph_edge_composite[end]) + for subnetwork_neighbor_node_id in subnetwork_neighbors_side_2 + if subnetwork_neighbor_node_id in keys(node_id_mapping) + allocgraph_node_id_2 = node_id_mapping[subnetwork_neighbor_node_id][1] + # Make sure this allocgraph node is distinct from the other one + if allocgraph_node_id_2 ≠ allocgraph_node_id_1 + push!(allocgraph_edge_composite, subnetwork_neighbor_node_id) + break + end + end + end + + # No connection to allocgraph node found on this side, so edge is discarded + if isnothing(allocgraph_node_id_2) + continue + end + + # Find capacity of this composite allocgraph edge + positive_flow = true + negative_flow = true + allocgraph_edge_capacity = Inf + # The start and end subnetwork nodes of the composite allocgraph + # edge are now nodes that have an equivalent in the allocgraph graph, + # these do not constrain the composite edge capacity + for allocgraph_edge_composite_triplet in + IterTools.partition(allocgraph_edge_composite, 3, 1) + node_type = lookup[allocgraph_edge_composite_triplet[2]] + node = getfield(p, node_type) + + # Find flow constraints + if is_flow_constraining(node) + problem_node_idx = + Ribasim.findsorted(node.node_id, allocgraph_edge_composite_triplet[2]) + allocgraph_edge_capacity = + min(allocgraph_edge_capacity, node.max_flow_rate[problem_node_idx]) + end + + # Find flow direction constraints + if is_flow_direction_constraining(node) + subnetwork_inneighbor_node_id = + only(inneighbors(graph_flow, allocgraph_edge_composite_triplet[2])) + + if subnetwork_inneighbor_node_id == allocgraph_edge_composite_triplet[1] + negative_flow = false + elseif subnetwork_inneighbor_node_id == allocgraph_edge_composite_triplet[3] + positive_flow = false + end + end + end + + # Add composite allocgraph edge(s) + if positive_flow + add_edge!(graph_allocation, allocgraph_node_id_1, allocgraph_node_id_2) + capacity[allocgraph_node_id_1, allocgraph_node_id_2] = allocgraph_edge_capacity + end + + if negative_flow + add_edge!(graph_allocation, allocgraph_node_id_2, allocgraph_node_id_1) + capacity[allocgraph_node_id_2, allocgraph_node_id_1] = allocgraph_edge_capacity + end + end + return capacity +end + +""" +The source nodes must only have one outneighbor. +""" +function valid_sources( + graph_allocation::DiGraph{Int}, + node_id_mapping::Dict{Int, Tuple{Int, Symbol}}, +)::Bool + errors = false + + for (allocgraph_node_id, allocgraph_node_type) in values(node_id_mapping) + if allocgraph_node_type == :source + if !( + (length(inneighbors(graph_allocation, allocgraph_node_id)) == 0) && + (length(outneighbors(graph_allocation, allocgraph_node_id)) == 1) + ) + @error "Sources nodes in the max flow graph must have no inneighbors and 1 outneighbor." + errors = true + end + end + end + return !errors +end + +""" +Remove user return flow edges that are upstream of the user itself, and collect the IDs +of the allocation graph node IDs of the users that do not have this problem. +""" +function avoid_using_own_returnflow!( + graph_allocation::DiGraph{Int}, + allocgraph_node_ids_user::Vector{Int}, + node_id_mapping_inverse::Dict{Int, Tuple{Int, Symbol}}, +)::Vector{Int} + allocgraph_node_ids_user_with_returnflow = Int[] + for allocgraph_node_id_user in allocgraph_node_ids_user + allocgraph_node_id_return_flow = + only(outneighbors(graph_allocation, allocgraph_node_id_user)) + if path_exists_in_graph( + graph_allocation, + allocgraph_node_id_return_flow, + allocgraph_node_id_user, + ) + rem_edge!( + graph_allocation, + allocgraph_node_id_user, + allocgraph_node_id_return_flow, + ) + @debug "The outflow of user #$(node_id_mapping_inverse[allocgraph_node_id_user][1]) is upstream of the user itself and thus ignored in allocation solves." + else + push!(allocgraph_node_ids_uer_with_returnflow, allocgraph_node_id_user) + end + end + return allocgraph_node_ids_user_with_returnflow +end + +""" +Build the graph used for the allocation problem. +""" +function allocation_graph( + p::Parameters, + subnetwork_node_ids::Vector{Int}, + source_edge_ids::Vector{Int}, +) + # Get the subnetwork and allocation node correspondence + node_id_mapping, source_edge_mapping = + get_node_id_mapping(p, subnetwork_node_ids, source_edge_ids) + + # Invert the node id mapping to easily translate from allocgraph nodes to subnetwork nodes + node_id_mapping_inverse = Dict{Int, Tuple{Int, Symbol}}() + for (subnetwork_node_id, (allocgraph_node_id, node_type)) in node_id_mapping + node_id_mapping_inverse[allocgraph_node_id] = (subnetwork_node_id, node_type) + end + + # Initialize the allocation graph + graph_allocation = DiGraph(length(node_id_mapping)) + + # Find the edges in the allocation graph + allocgraph_edges_composite, capacity = find_allocation_graph_edges!( + graph_allocation, + node_id_mapping, + p, + subnetwork_node_ids, + ) + + # Process the edges in the allocation graph + process_allocation_graph_edges!( + graph_allocation, + capacity, + allocgraph_edges_composite, + node_id_mapping, + p, + ) + + if !valid_sources(graph_allocation, node_id_mapping) + error("Errors in sources in allocation graph.") + end + + allocgraph_node_ids_user = [ + allocgraph_node_id for + (allocgraph_node_id, node_type) in values(node_id_mapping) if node_type == :user + ] + + allocgraph_node_ids_user_with_returnflow = avoid_using_own_returnflow!( + graph_allocation, + allocgraph_node_ids_user, + node_id_mapping_inverse, + ) + + # Used for updating user demand and source flow constraints + allocgraph_edges = collect(edges(graph_allocation)) + allocgraph_edge_ids_user_demand = Int[] + for (i, allocgraph_edge) in enumerate(allocgraph_edges) + allocgraph_node_type_dst = node_id_mapping_inverse[allocgraph_edge.dst][2] + if allocgraph_node_type_dst == :user + push!(allocgraph_edge_ids_user_demand, i) + end + end + + return graph_allocation, + capacity, + node_id_mapping, + node_id_mapping_inverse, + source_edge_mapping, + allocgraph_node_ids_user, + allocgraph_node_ids_user_with_returnflow, + allocgraph_edges, + allocgraph_edge_ids_user_demand +end + +""" +Add the flow variables F to the allocation problem. +The variable indices are the allocation graph edge IDs. +Non-negativivity constraints are also immediately added to the flow variables. +""" +function add_variables_flow!( + problem::JuMPModel, + allocgraph_edges::Vector{Edge{Int}}, +)::Nothing + n_flows = length(allocgraph_edges) + problem[:F] = @variable(problem, F[1:n_flows] >= 0.0) + return nothing +end + +""" +Add the user allocation variables A_user_{i} to the allocation problem. +The variable name indices i are the allocation graph user node IDs, +The variable indices are the priorities. +""" +function add_variables_allocation_user!( + problem::JuMPModel, + user::User, + allocgraph_edges::Vector{Edge{Int}}, + allocgraph_edge_ids_user_demand::Vector{Int}, +)::Nothing + for allocgraph_edge_id_user_demand in allocgraph_edge_ids_user_demand + allocgraph_node_id_user = allocgraph_edges[allocgraph_edge_id_user_demand].dst + base_name = "A_user_$allocgraph_node_id_user" + problem[Symbol(base_name)] = + @variable(problem, [1:length(user.priorities)], base_name = base_name) + end + return nothing +end + +""" +Add the basin allocation variables A_basin to the allocation problem. +The variable indices are the allocation graph basin node IDs. +Non-negativivity constraints are also immediately added to the basin allocation variables. +""" +function add_variables_allocation_basin!( + problem::JuMPModel, + node_id_mapping::Dict{Int, Tuple{Int, Symbol}}, + allocgraph_node_ids_basin::Vector{Int}, +)::Nothing + @variable(problem, A_basin[i = allocgraph_node_ids_basin] >= 0.0) + return nothing +end + +""" +Add the user allocation constraints to the allocation problem: +- The sum of the allocations to a user is equal to the flow to that user; +- The allocations to the users are non-negative; +- The allocations to the users are bounded from above by the user demands + (these are set before each allocation solve). + +The demand constrains have name demand_user_{i} where the i are the allocation graph +user node IDs and the constraint indices are the priorities. + +Constraints: +sum(allocations to user of all priorities) = flow to user +allocation to user at priority >= 0 +allocation to user at priority <= demand from user at priority +""" +function add_constraints_user_allocation!( + problem::JuMPModel, + user::User, + allocgraph_edges::Vector{Edge{Int}}, + allocgraph_edge_ids_user_demand::Vector{Int}, +)::Nothing + F = problem[:F] + for allocgraph_edge_id_user_demand in allocgraph_edge_ids_user_demand + allocgraph_node_id_user = allocgraph_edges[allocgraph_edge_id_user_demand].dst + base_name = "A_user_$allocgraph_node_id_user" + A_user = problem[Symbol(base_name)] + # Sum of allocations to user is total flow to user + @constraint( + problem, + sum(A_user) == F[allocgraph_edge_id_user_demand], + base_name = "allocation_sum[$allocgraph_node_id_user]" + ) + # Allocation flows are non-negative + @constraint(problem, [p = 1:length(user.priorities)], A_user[p] >= 0) + # Allocation flows are bounded from above by demands + base_name = "demand_user_$allocgraph_node_id_user" + problem[Symbol(base_name)] = @constraint( + problem, + [p = 1:length(user.priorities)], + A_user[p] <= 0, + base_name = base_name + ) + end + return nothing +end + +""" +Add the basin allocation constraints to the allocation problem; +the allocations to the basins are bounded from above by the basin demand +(these are set before each allocation solve). +The constraint indices are allocation graph basin node IDs. + +Constraint: +allocation to basin <= basin demand +""" +function add_constraints_basin_allocation!( + problem::JuMPModel, + allocgraph_node_ids_basin::Vector{Int}, +)::Nothing + A_basin = problem[:A_basin] + problem[:basin_allocation] = @constraint( + problem, + [i = allocgraph_node_ids_basin], + A_basin[i] <= 0.0, + base_name = "basin_allocation" + ) + return nothing +end + +""" +Add the flow capacity constraints to the allocation problem. +Only finite capacities get a constraint. +The constraint indices are the allocation graph edge IDs. + +Constraint: +flow over edge <= edge capacity +""" +function add_constraints_capacity!( + problem::JuMPModel, + capacity::SparseMatrixCSC{Float64, Int}, + allocgraph_edges::Vector{Edge{Int}}, +)::Nothing + F = problem[:F] + allocgraph_edge_ids_finite_capacity = Int[] + for (i, allocgraph_edge) in enumerate(allocgraph_edges) + if !isinf(capacity[allocgraph_edge.src, allocgraph_edge.dst]) + push!(allocgraph_edge_ids_finite_capacity, i) + end + end + problem[:capacity] = @constraint( + problem, + [i = allocgraph_edge_ids_finite_capacity], + F[i] <= capacity[allocgraph_edges[i].src, allocgraph_edges[i].dst], + base_name = "capacity" + ) + return nothing +end + +""" +Add the source constraints to the allocation problem. +The actual threshold values will be set before each allocation solve. +The constraint indices are the allocation graph source node IDs. + +Constraint: +flow over source edge <= source flow in subnetwork +""" +function add_constraints_source!( + problem::JuMPModel, + source_edge_mapping::Dict{Int, Int}, + allocgraph_edges::Vector{Edge{Int}}, + graph_allocation::DiGraph{Int}, +)::Nothing + F = problem[:F] + problem[:source] = @constraint( + problem, + [i = keys(source_edge_mapping)], + F[findfirst( + ==(Edge(i, only(outneighbors(graph_allocation, i)))), + allocgraph_edges, + )] <= 0.0, + base_name = "source" + ) + return nothing +end + +""" +Add the flow conservation constraints to the allocation problem. +The constraint indices are allocgraph user node IDs. + +Constraint: +sum(flows out of node node) <= flows into node + flow from storage and vertical fluxes +""" +function add_constraints_flow_conservation!( + problem::JuMPModel, + allocgraph_node_ids_basin::Vector{Int}, + allocgraph_node_inedge_ids::Dict{Int, Vector{Int}}, + allocgraph_node_outedge_ids::Dict{Int, Vector{Int}}, +)::Nothing + F = problem[:F] + problem[:flow_conservation] = @constraint( + problem, + [i = allocgraph_node_ids_basin], + sum([ + F[allocgraph_edge_id] for allocgraph_edge_id in allocgraph_node_outedge_ids[i] + ]) <= sum([ + F[allocgraph_edge_id] for allocgraph_edge_id in allocgraph_node_inedge_ids[i] + ]), + base_name = "flow_conservation", + ) + return nothing +end + +""" +Add the user returnflow constraints to the allocation problem. +The constraint indices are allocation graph user node IDs. + +Constraint: +outflow from user = return factor * inflow to user +""" +function add_constraints_user_returnflow!( + problem::JuMPModel, + allocgraph_node_ids_user_with_returnflow::Vector{Int}, +)::Nothing + F = problem[:F] + + problem[:return_flow] = @constraint( + problem, + [i = allocgraph_node_ids_user_with_returnflow], + F[only(allocgraph_node_outedge_ids[i])] == + user.return_factor[findsorted(user.node_id, node_id_mapping_inverse[i][1])] * + F[only(allocgraph_node_inedge_ids[i])], + base_name = "return_flow", + ) + return nothing +end + +""" +Add the objective function to the allocation problem. +Objective function: linear combination of allocations to the basins and users, where + basin allocations get a weight of 1.0 and user allocations get a weight of 2^(-priority index). +""" +function add_objective_function!( + problem::JuMPModel, + user::User, + allocgraph_node_ids_user::Vector{Int}, +)::Nothing + A_basin = problem[:A_basin] + allocation_user_weights = 1 ./ (2 .^ (1:length(user.priorities))) + allocation_user_variables = [ + problem[Symbol("A_user_$allocgraph_node_id_user")] for + allocgraph_node_id_user in allocgraph_node_ids_user + ] + @objective( + problem, + Max, + sum(A_basin) + sum([ + sum(allocation_user_variable .* allocation_user_weights) for + allocation_user_variable in allocation_user_variables + ]) + ) + return nothing +end + +""" +Construct the allocation problem for the current subnetwork as a JuMP.jl model. +""" +function allocation_problem( + p::Parameters, + node_id_mapping::Dict{Int, Tuple{Int, Symbol}}, + allocgraph_node_ids_user::Vector{Int}, + allocgraph_node_ids_user_with_returnflow::Vector{Int}, + allocgraph_edges::Vector{Edge{Int}}, + allocgraph_edge_ids_user_demand::Vector{Int}, + source_edge_mapping::Dict{Int, Int}, + graph_allocation::DiGraph{Int}, + capacity::SparseMatrixCSC{Float64, Int}, +)::JuMPModel + (; user) = p + + allocgraph_node_ids_basin = sort([ + allocgraph_node_id for + (allocgraph_node_id, node_type) in values(node_id_mapping) if node_type == :basin + ]) + + allocgraph_node_inedge_ids, allocgraph_node_outedge_ids = + get_node_in_out_edges(graph_allocation) + + problem = JuMPModel(HiGHS.Optimizer) + + # Add variables to problem + add_variables_flow!(problem, allocgraph_edges) + add_variables_allocation_user!( + problem, + user, + allocgraph_edges, + allocgraph_edge_ids_user_demand, + ) + add_variables_allocation_basin!(problem, node_id_mapping, allocgraph_node_ids_basin) + + # Add constraints to problem + add_constraints_user_allocation!( + problem, + user, + allocgraph_edges, + allocgraph_edge_ids_user_demand, + ) + add_constraints_basin_allocation!(problem, allocgraph_node_ids_basin) + add_constraints_capacity!(problem, capacity, allocgraph_edges) + add_constraints_source!( + problem, + source_edge_mapping, + allocgraph_edges, + graph_allocation, + ) + add_constraints_flow_conservation!( + problem, + allocgraph_node_ids_basin, + allocgraph_node_inedge_ids, + allocgraph_node_outedge_ids, + ) + add_constraints_user_returnflow!(problem, allocgraph_node_ids_user_with_returnflow) + # TODO: The fractional flow constraints + + # Add objective to problem + add_objective_function!(problem, user, allocgraph_node_ids_user) + + return problem +end + +""" +Construct the JuMP.jl problem for allocation. + +Definitions +----------- +- 'subnetwork' is used to refer to the original Ribasim subnetwork; +- 'allocgraph' is used to refer to the allocation graph. + +Inputs +------ +p: Ribasim problem parameters +subnetwork_node_ids: the problem node IDs that are part of the allocation subnetwork +source_edge_ids:: The IDs of the edges in the subnetwork whose flow fill be taken as + a source in allocation +Δt_allocation: The timestep between successive allocation solves + +Outputs +------- +An AllocationModel object. + +""" +function AllocationModel( + p::Parameters, + subnetwork_node_ids::Vector{Int}, + source_edge_ids::Vector{Int}, + Δt_allocation::Float64, +)::AllocationModel + graph_allocation, + capacity, + node_id_mapping, + node_id_mapping_inverse, + source_edge_mapping, + allocgraph_node_ids_user, + allocgraph_node_ids_user_with_returnflow, + allocgraph_edges, + allocgraph_edge_ids_user_demand = + allocation_graph(p, subnetwork_node_ids, source_edge_ids) + + # The JuMP.jl allocation problem + problem = allocation_problem( + p, + node_id_mapping, + allocgraph_node_ids_user, + allocgraph_node_ids_user_with_returnflow, + allocgraph_edges, + allocgraph_edge_ids_user_demand, + source_edge_mapping, + graph_allocation, + capacity, + ) + + return AllocationModel( + subnetwork_node_ids, + node_id_mapping, + node_id_mapping_inverse, + source_edge_mapping, + graph_allocation, + capacity, + problem, + Δt_allocation, + ) +end + +""" +Update the allocation problem with model data at the current: +- Demands of the users +- Flows of the source edges +- Demands of the basins +""" +function set_model_state_in_allocation!( + allocation_model::AllocationModel, + p::Parameters, + t::Float64, +)::Nothing + (; problem, node_id_mapping, source_edge_mapping) = allocation_model + (; user, connectivity) = p + (; priorities, demand) = user + (; flow, edge_ids_flow_inv) = connectivity + + # It is assumed that the allocation procedure does not have to be differentiated. + flow = get_tmp(flow, 0) + + for (subnetwork_node_id, (allocgraph_node_id, allocgraph_node_type)) in node_id_mapping + if allocgraph_node_type == :user + node_idx = findsorted(user.node_id, subnetwork_node_id) + demand_node = demand[node_idx] + base_name = "demand_user_$allocgraph_node_id" + constraints_demand = problem[Symbol(base_name)] + for priority_idx in eachindex(priorities) + set_normalized_rhs( + constraints_demand[priority_idx], + demand_node[priority_idx](t), + ) + end + elseif allocgraph_node_type == :source + subnetwork_edge = source_edge_mapping[allocgraph_node_id] + subnetwork_node_ids = edge_ids_flow_inv[subnetwork_edge] + constraint_source = problem[:source][allocgraph_node_id] + set_normalized_rhs(constraint_source, flow[subnetwork_node_ids]) + elseif allocgraph_node_type == :basin + # 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. + elseif allocgraph_node_type == :junction + nothing + else + error("Got unsupported allocation graph node type $allocgraph_node_type.") + end + end + return nothing +end + +""" +Assign the allocations to the users as determined by the solution of the allocation problem. +""" +function assign_allocations!(allocation_model::AllocationModel, user::User)::Nothing + (; problem, node_id_mapping) = allocation_model + for (subnetwork_node_id, (allocgraph_node_id, allocgraph_node_type)) in node_id_mapping + if allocgraph_node_type == :user + user_idx = findsorted(user.node_id, subnetwork_node_id) + base_name = "A_user_$allocgraph_node_id" + user.allocated[user_idx] .= value.(problem[Symbol(base_name)]) + 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)::Nothing + # Update allocation problem with data from main model + set_model_state_in_allocation!(allocation_model, p, t) + + # Solve the allocation problem + optimize!(allocation_model.problem) + + # Assign the allocations to the users + assign_allocations!(allocation_model, p.user) +end diff --git a/core/src/create.jl b/core/src/create.jl index e1e8b921a..6e0c0aa01 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -230,6 +230,9 @@ function Connectivity(db::DB, config::Config, chunk_size::Int)::Connectivity flow = FixedSizeDiffCache(flow, chunk_size) end + # TODO: Create allocation models from input here + allocation_models = AllocationModel[] + return Connectivity( graph_flow, graph_control, @@ -239,6 +242,7 @@ function Connectivity(db::DB, config::Config, chunk_size::Int)::Connectivity edge_ids_control, edge_connection_types_flow, edge_connection_types_control, + allocation_models, ) end diff --git a/core/src/solve.jl b/core/src/solve.jl index dd9bceb58..00522ea54 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -4,6 +4,30 @@ const ScalarInterpolation = const VectorInterpolation = LinearInterpolation{Vector{Vector{Float64}}, Vector{Float64}, true, Vector{Float64}} +""" +Store information for a subnetwork used for allocation. + +node_id: All the IDs of the nodes that are in this subnetwork +node_id_mapping: Mapping Dictionary; model_node_id => AG_node_id where such a correspondence exists + (all AG node ids are in the values) +node_id_mapping_inverse: The inverse of node_id_mapping, Dictionary; AG node ID => model node ID +Source edge mapping: AG source node ID => subnetwork source edge ID +graph_allocation: The graph used for the allocation problems +capacity: The capacity per edge of the allocation graph, 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 +""" +struct AllocationModel + node_id::Vector{Int} + node_id_mapping::Dict{Int, Tuple{Int, Symbol}} + node_id_mapping_inverse::Dict{Int, Tuple{Int, Symbol}} + source_edge_mapping::Dict{Int, Int} + graph_allocation::DiGraph{Int} + capacity::SparseMatrixCSC{Float64, Int} + problem::JuMPModel + Δt_allocation::Float64 +end + """ Store the connectivity information @@ -27,6 +51,7 @@ struct Connectivity{T} edge_ids_control::Dictionary{Tuple{Int, Int}, Int} edge_connection_type_flow::Dictionary{Int, Tuple{Symbol, Symbol}} edge_connection_type_control::Dictionary{Int, Tuple{Symbol, Symbol}} + allocation_models::Vector{AllocationModel} function Connectivity( graph_flow, graph_control, @@ -36,6 +61,7 @@ struct Connectivity{T} edge_ids_control, edge_connection_types_flow, edge_connection_types_control, + subnetwork, ) where {T} invalid_networks = Vector{String}() @@ -57,6 +83,7 @@ struct Connectivity{T} edge_ids_control, edge_connection_types_flow, edge_connection_types_control, + subnetwork, ) else invalid_networks = join(invalid_networks, ", ") diff --git a/core/src/utils.jl b/core/src/utils.jl index 70f9617eb..282007a1a 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -918,3 +918,48 @@ function reduction_factor(x::T, threshold::Real)::T where {T <: Real} one(T) end end + +"""Whether the given node node is flow constraining by having a maximum flow rate.""" +is_flow_constraining(node::AbstractParameterNode) = hasfield(typeof(node), :max_flow_rate) + +"""Whether the given node is flow direction constraining (only in direction of edges).""" +is_flow_direction_constraining(node::AbstractParameterNode) = + (nameof(typeof(node)) ∈ [:Pump, :Outlet, :TabulatedRatingCurve]) + +"""Find out whether a path exists between a start node and end node in the given graph.""" +function path_exists_in_graph(graph::DiGraph, start_node_id::Int, end_node_id::Int)::Bool + node_ids_visited = Set{Int}() + stack = [start_node_id] + + while !isempty(stack) + current_node_id = pop!(stack) + if current_node_id == end_node_id + return true + end + if !(current_node_id in node_ids_visited) + push!(node_ids_visited, current_node_id) + for outneighbor_node_id in outneighbors(graph, current_node_id) + push!(stack, outneighbor_node_id) + end + end + end + return false +end + +""" +Get two dictionaries, where: +- The first one gives the IDs of the inedges for each node ID in the graph +- The second one gives the IDs of the outedges for each node ID in the graph +""" +function get_node_in_out_edges( + graph::DiGraph{Int}, +)::Tuple{Dict{Int, Vector{Int}}, Dict{Int, Vector{Int}}} + n_nodes = nv(graph) + node_inedge_ids = Dict(i => Int[] for i in 1:n_nodes) + node_outedge_ids = Dict(i => Int[] for i in 1:n_nodes) + for (i, edge) in enumerate(edges(graph)) + push!(node_inedge_ids[edge.dst], i) + push!(node_outedge_ids[edge.src], i) + end + return node_inedge_ids, node_outedge_ids +end diff --git a/core/test/allocation.jl b/core/test/allocation.jl new file mode 100644 index 000000000..c92dd5dfa --- /dev/null +++ b/core/test/allocation.jl @@ -0,0 +1,34 @@ +import Ribasim +using SQLite +using JuMP: value + +@testset "Allocation solve" begin + toml_path = normpath(@__DIR__, "../../generated_testmodels/subnetwork/subnetwork.toml") + @test ispath(toml_path) + cfg = Ribasim.Config(toml_path) + gpkg_path = Ribasim.input_path(cfg, cfg.geopackage) + db = SQLite.DB(gpkg_path) + + p = Ribasim.Parameters(db, cfg) + close(db) + + # Inputs specific for this test model + subgraph_node_ids = unique(keys(p.lookup)) + source_edge_ids = [p.connectivity.edge_ids_flow[(1, 2)]] + flow = Ribasim.get_tmp(p.connectivity.flow, 0) + flow[1, 2] = 4.5 # Source flow + Δt_allocation = 24.0 * 60^2 + t = 0.0 + + allocation_model = + Ribasim.AllocationModel(p, subgraph_node_ids, source_edge_ids, Δt_allocation) + Ribasim.allocate!(p, allocation_model, t) + + F = value.(allocation_model.problem[:F]) + @test F ≈ [3.0, 3.0, 0.5, 3.5, 1.0, 4.5] + + allocated = p.user.allocated + @test allocated[1] ≈ [0.0, 1.0, 0.0] + @test allocated[2] ≈ [0.0, 0.5, 0.0] + @test allocated[3] ≈ [3.0, 0.0, 0.0] +end diff --git a/core/test/runtests.jl b/core/test/runtests.jl index 8a73f26bc..6d6c293ff 100644 --- a/core/test/runtests.jl +++ b/core/test/runtests.jl @@ -12,6 +12,7 @@ using SafeTestsets: @safetestset @safetestset "Basic Model Interface" include("bmi.jl") @safetestset "Utility functions" include("utils.jl") @safetestset "Control" include("control.jl") + @safetestset "Allocation" include("allocation.jl") @safetestset "Time" include("time.jl") @safetestset "Docs" include("docs.jl") @safetestset "Command Line Interface" include("cli.jl") diff --git a/core/test/time.jl b/core/test/time.jl index 6724ea808..da8a996ea 100644 --- a/core/test/time.jl +++ b/core/test/time.jl @@ -45,7 +45,6 @@ end t_end = Ribasim.seconds_since(cfg.endtime, cfg.starttime) - # demand[user_idx][priority](t) @test demand[1][2](0.5 * t_end) ≈ 1.0 @test demand[2][1](0.0) ≈ 0.0 @test demand[2][1](t_end) ≈ 0.0 diff --git a/docs/_quarto.yml b/docs/_quarto.yml index 9f7b2891a..5025f5a15 100644 --- a/docs/_quarto.yml +++ b/docs/_quarto.yml @@ -30,6 +30,7 @@ website: - core/index.qmd - core/usage.qmd - core/equations.qmd + - core/allocation.qmd - core/numerics.qmd - build/index.md - title: "Python tooling" diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd new file mode 100644 index 000000000..a8706f8d6 --- /dev/null +++ b/docs/core/allocation.qmd @@ -0,0 +1,223 @@ +--- +title: "Allocation" +--- + +Allocation is the process of assigning an allocated abstraction flow rate to user nodes in the model based on information about sources, user demands over various priorities, constraints introduced by nodes, local water availability and graph topology. The allocation procedure implemented in Ribasim is heavily inspired by the [maximum flow problem](https://en.wikipedia.org/wiki/Maximum_flow_problem). + +The allocation problem is solved per subnetwork of the Ribasim model. The subnetwork is used to formulate an optimization problem with the [JuMP](https://jump.dev/JuMP.jl/stable/) package, which is solved using the [HiGHS solver](https://highs.dev/). See also the example of solving the maximum flow problem with `JuMP.jl` [here](https://jump.dev/JuMP.jl/stable/tutorials/linear/network_flows/#The-max-flow-problem). + +# The allocation problem + +The following data of the parameters and state of a Ribasim model are relevant for the allocation problem. + +## Allocation problem input + +### The subnetwork + +The allocation problem is solved per subgraph, where a subgraph is given by a subset $S \subset V$ of node ids. Different subgraphs are disjoint from eachother. + +### Source flows + +Sources are indicated by a set of edges in the subnetwork +$$ +E_S^\text{source} \subset \left(S \times S\right) \cap E. +$$ +That is, if $(i,j) \in E_S^\text{source}$, then $Q_{ij}$ is treated as a source flow in the allocation problem. + +### User demands + +The subnetwork contains a subset of user nodes $U_S \subset S$, who all have time varying demands over various priorities $p$: +$$ + d^p_i(t), \quad i \in U_S, p = 1,2,\ldots, p_{\max}. +$$ + +:::{.callout-note} +On this page we assume that the priorities are given by all integers from $1$ to some $p_{\max} \in \mathbb{N}$. However, in the Ribasim input this is not a requirement; some of these in between priority values can be missing, only the ordering of the given priorities is taken into account. +::: + +### 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 +$$ + \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}$ +$$ + u_i(t)-S(l_{\min,i}), \quad \forall i \in B_S. +$$ +Note that this value can be negative, which we interpret as a demand from the basin. + +### Flow magnitude and direction constraints + +Nodes in the Ribasim model that have a `max_flow_rate`, i.e. pumps and outlets, put a constraint on the flow through that node. Some nodes only allow flow in one direction, like pumps, outlets and tabulated rating curves. + +### Fractional flows and user return flows + +Both fractional flow nodes and user nodes dictate proportional relationships between flows over edges in the subnetwork. Users have a return factor $r_i, i \in U_S$. + +## The allocation optimization problem + +### The allocation graph + +A new graph is created from the subnetwork, which we call the allocation graph. To indicate the difference between subnetwork data and allocation graph data, the allocation graph data is denoted with a hat. The allocation graph consists of: + +- Nodes $\hat{V_S}$, where each basin, source and user in the subnetwork get a node in the allocation graph. Also nodes that have fractional flow outneighbors get a node in the allocation graph. The implementation makes heavy use of the node id mapping $m_S : i \mapsto \hat{i}$ to translate from subnetwork node IDs to allocation graph node IDs. Unless specified otherwise, we assume this relationship between index symbols that appear both with and without a hat. +- Edges $\hat{E_S}$, where the edges in the allocation graph are given by one or more edges in the subnetwork, where those edges connect nodes in the subnetwork that have an equivalent in the allocation graph. The direction of the edges in the allocation graph is given by the direction constraints in the subnetwork. + +For notational convenience, we use the notation +$$ + \begin{align} + \hat{V}^{\text{out}}_S(\hat{i}) = \left\{\hat{j} \in \hat{V}_S : (\hat{i},\hat{j}) \in \hat{E}_S\right\} \\ + \hat{V}^{\text{in}}_S(\hat{j}) = \left\{\hat{i} \in \hat{V}_S : (\hat{i},\hat{j}) \in \hat{E}_S\right\} + \end{align} +$$ +for the set of in-neighbors and out-neighbors of a node in the allocation graph respectively. + +### The allocation graph capacities + +The capacities of the edges of the allocation graph are collected in the sparse capacity matrix $\hat{C}_S \in \overline{\mathbb{R}}_{\ge 0}^{\hat{n}\times\hat{n}}$ where $\hat{n}$ is the number of nodes in the allocation graph. The capacities can be infinite. + +The capacities are determined in 3 different ways: + +- If an edge does not exist, i.e. $(\hat{i},\hat{j}) \notin \hat{E}$ for certain $1 \le \hat{i},\hat{j}\le \hat{n}$, then $(\hat{C}_S)_{\hat{i},\hat{j}} = 0$; +- The capacity of the edge $\hat{e} \in \hat{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. + +### The optimization variables + +There are several types of variables whose value has to be determined to solve the allocation problem: + +- The flows $F \in \mathbb{R}_{\ge 0}^{\hat{n}\times\hat{n}}$ over the edges in the allocation graph; +- The allocations to the users +$$ + A^\text{user}_{\hat{i},p} \ge 0, \quad \forall \hat{i} \in \hat{U}_S, \forall p \in \{1,2,\ldots, p_\max\}, +$$ +where $\hat{U}_S = m_S(U_S) \subset \hat{V}_S$ is the set of user node ids in the allocation graph; + +- The allocations to the basins +$$ + A^\text{basin}_{\hat{i}} \ge 0, \quad \hat{B}_S, +$$ +where $\hat{B} = m_S(B_S) \subset \hat{V}_S$ is the set of basin node ids in the allocation graph. + +### The optimization objective + +The goal of allocation is to maximize the flow to the users. However, basins can also demand water if their level is below the minimum abstraction level, which we give a higher priority than user demands. Therefore, we use the following optimization objective: +$$ + \max \sum_{\hat{i} \in \hat{B}_S} A^\text{basin}_{\hat{i}} + \sum_{\hat{i}\in \hat{U}_S} \sum_{p=1}^{p_\max} 2^{-p} A^\text{user}_{\hat{i},p}. +$$ {#eq-objective} +This is a linear combination of all allocations, where allocations to basins get a weight of $1$ and allocations to users get a weight of $2^{-p}$ where $p$ is the priority. + +### The optimization variable constraints + +- Source flows: for the source edges we have that +$$ +F_{\hat{i}\hat{j}} \le Q_{ij} \quad \forall (i,j) \in E_S^\text{source}. +$$ +Note that we do not require equality here; not all source flow has to be used. +- Flow conservation: For the basins in the allocation graph we have that +$$ + \sum_{\hat{j}=1}^{\hat{n}} F_{\hat{k}\hat{j}} \le \sum_{\hat{i}=1}^{\hat{n}} F_{\hat{i}\hat{k}} + \Phi_{\hat{k}}, \quad \forall\hat{k} \in \hat{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. +$\Phi_{\hat{k}}$ is the local water supply of the basins: +$$ + \Phi_{\hat{k}} = \max\left(\phi_k(t) + \frac{u_k(t)-S(l_{\min,k})}{\Delta t_\text{alloc}}, 0.0 \right). +$$ {#eq-basinsourceflow} +Here the first term denotes the vertical fluxes and the second term the flow that can be supplied by the water in the basin above its minimum level, where $\Delta t_\text{alloc}$ is the allocation solve timestep. + +- Capacity: the flows over the edges are positive and bounded by the edge capacity: +$$ + F_{\hat{i}\hat{j}} \le \left(\hat{C}_S\right)_{\hat{i}\hat{j}}, \quad \forall(\hat{i},\hat{j}) \in \hat{E}_S. +$$ {#eq-capacityconstraint} +- User outflow: The outflow of the user is dictated by the inflow and the return factor: +$$ +\begin{align} + F_{\hat{i}\hat{k}} = r_k \cdot F_{\hat{k}\hat{j}}\\ + \quad \forall\hat{k} \in \hat{U}_s, \\ + \hat{V}^{\text{in}}_S(\hat{k}) = \{\hat{i}\},\\ + \hat{V}^{\text{out}}_S(\hat{k}) = \{\hat{j}\}. +\end{align} +$$ {#eq-returnflowconstraint} +Here we use that each user node in the allocation graph has an unique in-edge and out-edge. +- User allocation: The flow over the edge to the user is equal to the sum of the allocations to the user: +$$ + F_{\hat{i}\hat{k}} = \sum_{p=1}^{p_\max} A^\text{user}_{\hat{k},p}, \quad \forall \hat{k} \in \hat{U}_S, \hat{V}^{\text{out}}_s(\hat{k}) = \{\hat{i}\}. +$$ {#eq-userallocationconstraint} +Here we use that each user has an unique out-edge. +- User demand: what is allocated to the user is bounded above by the user demand: +$$ + A_{\hat{i},p}^\text{user} \leq d_i^p(t) \quad \forall\hat{i} \in \hat{U}_S, \; p = 1,\ldots,p_\max. +$$ +- Basin allocation: If the flow supplied by a basin (as determined in @eq-basinsourceflow) is negative, it is a demand: +$$ +A^\text{basin}_{\hat{i}} = \max\left(-\left(\phi_i(t) + \frac{u_i(t)-S(l_{\min,k})}{\Delta t_\text{alloc}}\right), 0.0 \right), \quad \forall \hat{i} \in \hat{B}_S. +$$ {#eq-basinallocationconstraint} +- Fractinal flow: Let $\hat{L}_S \subset \hat{V}_S$ be the set of nodes in the max flow graph with fractional flow outneighbors, and $f_j$ the flow fraction associated with fractional flow node $j \in V_S$. Then +$$ +\begin{align} + F_{\hat{i}\hat{j}} = f_j \sum_{k\in \hat{V}^\text{in}_S(\hat{i})} F_{\hat{k}\hat{i}} \\ + \forall \hat{i} \in \hat{L}_S, \\ + \hat{j} \in \hat{V}_S^\text{out}(\hat{i}). +\end{align} +$$ {#eq-fractionalflowconstraint} + +- Flow sign: Furthermore there are the non-negativity constraints for the flows and allocations, see [The optimization variables](allocation.qmd#the-optimization-variables). + +## Final notes on the allocation problem + +### Users using their own return flow + +If not explicitly avoided, users can use their own return flow in this allocation problem formulation. +Therefore, return flow of users is only taken into account by allocation if that return flow is downstream of the user where it comes from. That is, if there is no path in the directed allocation graph from the user outflow node back to the user. + +# Solving the allocation problem + +The text below shows a representation generated by `JuMP.jl` of an optimization as described in the previous section. + +``` +Max A_basin[2] + A_basin[3] + A_basin[5] + 0.5 A_user_4[1] + 0.25 A_user_4[2] + 0.125 A_user_4[3] + 0.5 A_user_6[1] + 0.25 A_user_6[2] + 0.125 A_user_6[3] + 0.5 A_user_1[1] + 0.25 A_user_1[2] + 0.125 A_user_1[3] +Subject to + allocation_sum[1] : -F[1] + A_user_1[1] + A_user_1[2] + A_user_1[3] == 0 + allocation_sum[4] : -F[3] + A_user_4[1] + A_user_4[2] + A_user_4[3] == 0 + allocation_sum[6] : -F[5] + A_user_6[1] + A_user_6[2] + A_user_6[3] == 0 + A_user_1[1] >= 0 + A_user_1[2] >= 0 + A_user_1[3] >= 0 + A_user_4[1] >= 0 + A_user_4[2] >= 0 + A_user_4[3] >= 0 + A_user_6[1] >= 0 + A_user_6[2] >= 0 + A_user_6[3] >= 0 + demand_user_1[1] : A_user_1[1] <= 4 + demand_user_1[2] : A_user_1[2] <= 0 + demand_user_1[3] : A_user_1[3] <= 0 + demand_user_4[1] : A_user_4[1] <= 0 + demand_user_4[2] : A_user_4[2] <= 2 + demand_user_4[3] : A_user_4[3] <= 0 + demand_user_6[1] : A_user_6[1] <= 0 + demand_user_6[2] : A_user_6[2] <= 1 + demand_user_6[3] : A_user_6[3] <= 0 + basin_allocation[2] : A_basin[2] <= 0 + basin_allocation[3] : A_basin[3] <= 0 + basin_allocation[5] : A_basin[5] <= 0 + capacity[2] : F[2] <= 3 + capacity[4] : F[4] <= 4 + source[7] : F[6] <= 4.5 + flow_conservation[2] : F[1] - F[2] <= 0 + flow_conservation[3] : F[2] + F[3] - F[4] <= 0 + flow_conservation[5] : F[4] + F[5] - F[6] <= 0 + F[1] >= 0 + F[2] >= 0 + F[3] >= 0 + F[4] >= 0 + F[5] >= 0 + F[6] >= 0 + A_basin[2] >= 0 + A_basin[3] >= 0 + A_basin[5] >= 0 +``` + +A more detailed explanation of this will follow in the future. diff --git a/docs/core/index.qmd b/docs/core/index.qmd index 8f06cd9e0..407dd2bf7 100644 --- a/docs/core/index.qmd +++ b/docs/core/index.qmd @@ -5,7 +5,7 @@ title: "Julia core" With the term "core", we mean the computational engine of Ribasim. As detailed in the [usage](usage.qmd) documentation, it is generally used as a command line tool. -The theory is described on the [equations](equations.qmd) page, and more in-depth numerical considerations are described on the [numerical considerations](numerics.qmd) page. +The theory is described on the [equations](equations.qmd) page, and more in-depth numerical considerations are described on the [numerical considerations](numerics.qmd) page. As allocation is a large and self-contained part of the Ribasim core, it is described on the separate [allocation](allocation.qmd) page. The core is implemented in the [Julia programming language](https://julialang.org/), and can be found in the [Ribasim repository](https://github.com/Deltares/Ribasim) under the diff --git a/python/ribasim/ribasim/node_types/user.py b/python/ribasim/ribasim/node_types/user.py index 258b8b910..2dda210b5 100644 --- a/python/ribasim/ribasim/node_types/user.py +++ b/python/ribasim/ribasim/node_types/user.py @@ -44,7 +44,9 @@ class Config: def sort(self): if self.static is not None: - self.static.sort_values("node_id", ignore_index=True, inplace=True) + self.static.sort_values( + ["node_id", "priority"], ignore_index=True, inplace=True + ) if self.time is not None: self.time.sort_values( ["node_id", "priority", "time"], ignore_index=True, inplace=True diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index 99b71a67f..028a1d947 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -1,6 +1,10 @@ __version__ = "0.2.0" -from ribasim_testmodels.allocation import subnetwork_model, user_model +from ribasim_testmodels.allocation import ( + looped_subnetwork_model, + subnetwork_model, + user_model, +) from ribasim_testmodels.backwater import backwater_model from ribasim_testmodels.basic import ( basic_model, @@ -65,4 +69,5 @@ "outlet_model", "user_model", "subnetwork_model", + "looped_subnetwork_model", ] diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index dffbd673d..2c9b6786e 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -207,15 +207,15 @@ def subnetwork_model(): } ) - # Setup the flow boundary: - flow_boundary = ribasim.FlowBoundary( - static=pd.DataFrame(data={"node_id": [1], "flow_rate": [1e-3]}) - ) - state = pd.DataFrame(data={"node_id": [2, 6, 8], "level": 1.0}) basin = ribasim.Basin(profile=profile, static=static, state=state) + # Setup the flow boundary: + flow_boundary = ribasim.FlowBoundary( + static=pd.DataFrame(data={"node_id": [1], "flow_rate": [4.5]}) + ) + # Setup the users: user = ribasim.User( static=pd.DataFrame( @@ -250,13 +250,16 @@ def subnetwork_model(): data={ "node_id": [5], "flow_rate": [4.0], + "max_flow_rate": [4.0], } ) ) # Setup the outlets: outlet = ribasim.Outlet( - static=pd.DataFrame(data={"node_id": [3, 7, 13], "flow_rate": 3.0}) + static=pd.DataFrame( + data={"node_id": [3, 7, 13], "flow_rate": 3.0, "max_flow_rate": 3.0} + ) ) # Setup the terminal: @@ -283,3 +286,252 @@ def subnetwork_model(): ) return model + + +def looped_subnetwork_model(): + # Setup the nodes: + xy = np.array( + [ + (0.0, 0.0), # 1: User + (0.0, 1.0), # 2: Basin + (-1.0, 1.0), # 3: Outlet + (-2.0, 1.0), # 4: Terminal + (2.0, 1.0), # 5: FlowBoundary + (0.0, 2.0), # 6: Pump + (-2.0, 3.0), # 7: Basin + (-1.0, 3.0), # 8: Outlet + (0.0, 3.0), # 9: Basin + (1.0, 3.0), # 10: Outlet + (2.0, 3.0), # 11: Basin + (-2.0, 4.0), # 12: User + (0.0, 4.0), # 13: TabulatedRatingCurve + (2.0, 4.0), # 14: TabulatedRatingCurve + (0.0, 5.0), # 15: Basin + (1.0, 5.0), # 16: Pump + (2.0, 5.0), # 17: Basin + (-1.0, 6.0), # 18: User + (0.0, 6.0), # 19: TabulatedRatingCurve + (2.0, 6.0), # 20: User + (0.0, 7.0), # 21: Basin + (0.0, 8.0), # 22: Outlet + (0.0, 9.0), # 23: Terminal + (3.0, 3.0), # 24: User + ] + ) + node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1]) + + node_type = [ + "User", + "Basin", + "Outlet", + "Terminal", + "FlowBoundary", + "Pump", + "Basin", + "Outlet", + "Basin", + "Outlet", + "Basin", + "User", + "TabulatedRatingCurve", + "TabulatedRatingCurve", + "Basin", + "Pump", + "Basin", + "User", + "TabulatedRatingCurve", + "User", + "Basin", + "Outlet", + "Terminal", + "User", + ] + + # Make sure the feature id starts at 1: explicitly give an index. + node = ribasim.Node( + static=gpd.GeoDataFrame( + data={"type": node_type}, + index=pd.Index(np.arange(len(xy)) + 1, name="fid"), + geometry=node_xy, + crs="EPSG:28992", + ) + ) + + # Setup the edges: + from_id = np.array( + [ + 5, + 2, + 3, + 2, + 2, + 6, + 9, + 8, + 7, + 9, + 13, + 15, + 16, + 17, + 15, + 19, + 15, + 18, + 21, + 22, + 9, + 10, + 11, + 14, + 1, + 12, + 20, + 11, + 24, + ], + dtype=np.int64, + ) + to_id = np.array( + [ + 2, + 3, + 4, + 1, + 6, + 9, + 8, + 7, + 12, + 13, + 15, + 16, + 17, + 20, + 19, + 21, + 18, + 21, + 22, + 23, + 10, + 11, + 14, + 17, + 2, + 7, + 17, + 24, + 11, + ], + dtype=np.int64, + ) + lines = ribasim.utils.geometry_from_connectivity(node, from_id, to_id) + edge = ribasim.Edge( + static=gpd.GeoDataFrame( + data={ + "from_node_id": from_id, + "to_node_id": to_id, + "edge_type": len(from_id) * ["flow"], + }, + geometry=lines, + crs="EPSG:28992", + ) + ) + + # Setup the basins: + profile = pd.DataFrame( + data={ + "node_id": [2, 2, 7, 7, 9, 9, 11, 11, 15, 15, 17, 17, 21, 21], + "area": 1000.0, + "level": 7 * [0.0, 1.0], + } + ) + + static = pd.DataFrame( + data={ + "node_id": [2, 7, 11, 9, 15, 17, 21], + "drainage": 0.0, + "potential_evaporation": 0.0, + "infiltration": 0.0, + "precipitation": 0.0, + "urban_runoff": 0.0, + } + ) + + state = pd.DataFrame(data={"node_id": [2, 7, 9, 15, 17, 21], "level": 1.0}) + + basin = ribasim.Basin(profile=profile, static=static, state=state) + + # Setup the flow boundary: + flow_boundary = ribasim.FlowBoundary( + static=pd.DataFrame(data={"node_id": [5], "flow_rate": [4.5]}) + ) + + # Setup the users: + user = ribasim.User( + static=pd.DataFrame( + data={ + "node_id": [1, 12, 18, 20, 24], + "demand": 1.0, + "return_factor": 0.9, + "min_level": 0.9, + "priority": [2, 1, 3, 3, 2], + } + ) + ) + + # Setup the pumps: + pump = ribasim.Pump( + static=pd.DataFrame( + data={ + "node_id": [6, 16], + "flow_rate": 4.0, + "max_flow_rate": 4.0, + } + ) + ) + + # Setup the outlets: + outlet = ribasim.Outlet( + static=pd.DataFrame( + data={"node_id": [3, 8, 10, 22], "flow_rate": 3.0, "max_flow_rate": 3.0} + ) + ) + + # Setup the tabulated rating curves + rating_curve = ribasim.TabulatedRatingCurve( + static=pd.DataFrame( + data={ + "node_id": [13, 13, 14, 14, 19, 19], + "level": 3 * [0.0, 1.0], + "discharge": 3 * [0.0, 2.0], + } + ) + ) + + # Setup the terminals: + terminal = ribasim.Terminal( + static=pd.DataFrame( + data={ + "node_id": [4, 23], + } + ) + ) + + model = ribasim.Model( + modelname="looped_subnetwork", + node=node, + edge=edge, + basin=basin, + flow_boundary=flow_boundary, + user=user, + pump=pump, + outlet=outlet, + tabulated_rating_curve=rating_curve, + terminal=terminal, + starttime="2020-01-01 00:00:00", + endtime="2021-01-01 00:00:00", + ) + + return model diff --git a/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py b/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py index f0942de6e..73551d2d1 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py +++ b/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py @@ -292,16 +292,3 @@ def dutch_waterways_model(): ) return model - - -if __name__ == "__main__": - import matplotlib.pyplot as plt - - model = dutch_waterways_model() - - model.plot() - - df_flow = model.flow_boundary.time.pivot_table(index="time", values=["flow_rate"]) - df_flow.plot() - - plt.show()