From 793fdd788cbd1a47348dfd886755bf37ee7478b7 Mon Sep 17 00:00:00 2001 From: Hofer-Julian <30049909+Hofer-Julian@users.noreply.github.com> Date: Fri, 2 Feb 2024 10:54:37 +0100 Subject: [PATCH 1/4] Move release script to Python and run it with pixi (#1029) Fixes #987 --------- Co-authored-by: Martijn Visser --- .../Ribasim_Ribasim_MakeGitHubRelease.xml | 25 +------ pixi.lock | 74 +++++++++++++++++++ pixi.toml | 3 + utils/github-release.py | 24 ++++++ 4 files changed, 104 insertions(+), 22 deletions(-) create mode 100644 utils/github-release.py diff --git a/.teamcity/Ribasim/buildTypes/Ribasim_Ribasim_MakeGitHubRelease.xml b/.teamcity/Ribasim/buildTypes/Ribasim_Ribasim_MakeGitHubRelease.xml index 3dd85bf1c..341866805 100644 --- a/.teamcity/Ribasim/buildTypes/Ribasim_Ribasim_MakeGitHubRelease.xml +++ b/.teamcity/Ribasim/buildTypes/Ribasim_Ribasim_MakeGitHubRelease.xml @@ -13,27 +13,9 @@ set -euxo pipefail . /usr/share/Modules/init/bash -module load github -# Get the name of the currently checked out tag -tag_name=$(git describe --tags --exact-match 2>/dev/null) - -# Check if a tag is checked out -if [ -n "$tag_name" ]; then - echo "Currently checked out tag: $tag_name" - - # Create a release using gh - gh release create "$tag_name" \ - --generate-notes \ - ribasim_cli_linux.zip \ - ribasim_cli_windows.zip \ - ribasim_qgis.zip \ - generated_testmodels.zip - - echo "Release created successfully." - -else - echo "No tag is currently checked out." -fi]]> +module load pixi +pixi run github-release +]]> @@ -130,4 +112,3 @@ fi]]> - diff --git a/pixi.lock b/pixi.lock index b076825ff..97d5014bb 100644 --- a/pixi.lock +++ b/pixi.lock @@ -9869,6 +9869,80 @@ package: license_family: BSD size: 82263 timestamp: 1594304231182 +- platform: linux-64 + name: gh + version: 2.43.1 + category: main + manager: conda + dependencies: [] + url: https://conda.anaconda.org/conda-forge/linux-64/gh-2.43.1-ha8f183a_0.conda + hash: + md5: adaa9131473d962cc1521d2169d59c60 + sha256: f7a7918096362a9f197c1b043a6e1d524fc8e077be6797972be910b4482c7693 + build: ha8f183a_0 + arch: x86_64 + subdir: linux-64 + build_number: 0 + license: Apache-2.0 + license_family: APACHE + size: 16858494 + timestamp: 1706752640334 +- platform: osx-64 + name: gh + version: 2.43.1 + category: main + manager: conda + dependencies: [] + url: https://conda.anaconda.org/conda-forge/osx-64/gh-2.43.1-h990441c_0.conda + hash: + md5: 620d1c108bc7f37c22102d181816eb19 + sha256: 8087704d81711ea3e9254d84cab3cf560ac2e79a668e2e2f9d875f1053707e70 + build: h990441c_0 + arch: x86_64 + subdir: osx-64 + build_number: 0 + constrains: + - __osx>=10.12 + license: Apache-2.0 + license_family: APACHE + size: 16650708 + timestamp: 1706752915101 +- platform: osx-arm64 + name: gh + version: 2.43.1 + category: main + manager: conda + dependencies: [] + url: https://conda.anaconda.org/conda-forge/osx-arm64/gh-2.43.1-h75b854d_0.conda + hash: + md5: 2fafc8df33219ebcc5a3eb69afef5f12 + sha256: 34783d293ab1db89469faa156e0c8e28fd213fa1d2414fdaaa008f7e4e31b45e + build: h75b854d_0 + arch: aarch64 + subdir: osx-arm64 + build_number: 0 + license: Apache-2.0 + license_family: APACHE + size: 16033884 + timestamp: 1706752975233 +- platform: win-64 + name: gh + version: 2.43.1 + category: main + manager: conda + dependencies: [] + url: https://conda.anaconda.org/conda-forge/win-64/gh-2.43.1-hd02998f_0.conda + hash: + md5: 229cbf4b096a468d30085cb9eead7a20 + sha256: 210b5dc8dcf91a9d528922b0ed9697a75a332f6c501dc95b806ee0f1ee6eda10 + build: hd02998f_0 + arch: x86_64 + subdir: win-64 + build_number: 0 + license: Apache-2.0 + license_family: APACHE + size: 16789533 + timestamp: 1706753295429 - platform: linux-64 name: giflib version: 5.2.1 diff --git a/pixi.toml b/pixi.toml index e29bfb276..f0f7583fe 100644 --- a/pixi.toml +++ b/pixi.toml @@ -122,11 +122,14 @@ install-qgis-plugins = { depends_on = [ mypy-ribasim-qgis = "mypy ribasim_qgis" # Run ribasim-model = "julia --project=core -e 'using Ribasim; Ribasim.main(ARGS)'" +# Release +github-release = "python utils/github-release.py" [dependencies] build = "*" datamodel-code-generator = "0.24.*" geopandas = "*" +gh = "*" hatchling = "*" juliaup = "*" jupyterlab = "*" diff --git a/utils/github-release.py b/utils/github-release.py new file mode 100644 index 000000000..c0d3da74e --- /dev/null +++ b/utils/github-release.py @@ -0,0 +1,24 @@ +import subprocess + +# Get the name of the currently checked out tag +tag_name = subprocess.check_output( + ["git", "describe", "--tags", "--exact-match"], text=True +).strip() + + +print(f"Currently checked out tag: {tag_name}") + +# Create a release using gh +subprocess.check_call( + [ + "gh", + "release", + "create", + tag_name, + "--generate-notes", + "ribasim_cli_linux.zip", + "ribasim_cli_windows.zip", + "ribasim_qgis.zip", + "generated_testmodels.zip", + ] +) From 2dc5d7478252976cf9569694475c83bf0c2b0467 Mon Sep 17 00:00:00 2001 From: Huite Date: Fri, 2 Feb 2024 10:58:49 +0100 Subject: [PATCH 2/4] Allow missing 'Basin / time' data (#1028) Fixes #1027 Allows for missing values in basin.time data. --------- Co-authored-by: Martijn Visser --- core/src/create.jl | 8 +- core/src/utils.jl | 6 +- core/src/validation.jl | 10 +-- core/test/run_models_test.jl | 36 +++++++++ docs/schema/BasinTime.schema.json | 62 +++++++++++---- python/ribasim/ribasim/models.py | 10 +-- .../ribasim_testmodels/__init__.py | 43 +++++------ .../ribasim_testmodels/bucket.py | 75 +++++++++++++++++++ 8 files changed, 196 insertions(+), 54 deletions(-) diff --git a/core/src/create.jl b/core/src/create.jl index 66e331e06..331a944f3 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -485,10 +485,10 @@ function Basin(db::DB, config::Config, chunk_sizes::Vector{Int})::Basin current_area = DiffCache(current_area, chunk_sizes) end - precipitation = fill(NaN, length(node_id)) - potential_evaporation = fill(NaN, length(node_id)) - drainage = fill(NaN, length(node_id)) - infiltration = fill(NaN, length(node_id)) + precipitation = zeros(length(node_id)) + potential_evaporation = zeros(length(node_id)) + drainage = zeros(length(node_id)) + infiltration = zeros(length(node_id)) table = (; precipitation, potential_evaporation, drainage, infiltration) area, level, storage = create_storage_tables(db, config) diff --git a/core/src/utils.jl b/core/src/utils.jl index 7d2bd8e64..27632f41b 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -625,12 +625,12 @@ end Update `table` at row index `i`, with the values of a given row. `table` must be a NamedTuple of vectors with all variables that must be loaded. The row must contain all the column names that are present in the table. -If a value is NaN, it is not set. +If a value is missing, it is not set. """ function set_table_row!(table::NamedTuple, row, i::Int)::NamedTuple for (symbol, vector) in pairs(table) val = getproperty(row, symbol) - if !isnan(val) + if !ismissing(val) vector[i] = val end end @@ -672,7 +672,7 @@ function set_current_value!( for (i, id) in enumerate(node_id) for (symbol, vector) in pairs(table) idx = findlast( - row -> row.node_id == id && !isnan(getproperty(row, symbol)), + row -> row.node_id == id && !ismissing(getproperty(row, symbol)), pre_table, ) if idx !== nothing diff --git a/core/src/validation.jl b/core/src/validation.jl index f74d03eb8..98d99d03c 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -185,11 +185,11 @@ end @version BasinTimeV1 begin node_id::Int time::DateTime - drainage::Float64 - potential_evaporation::Float64 - infiltration::Float64 - precipitation::Float64 - urban_runoff::Float64 + drainage::Union{Missing, Float64} + potential_evaporation::Union{Missing, Float64} + infiltration::Union{Missing, Float64} + precipitation::Union{Missing, Float64} + urban_runoff::Union{Missing, Float64} end @version BasinProfileV1 begin diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index 109c20dba..b2d992454 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -98,6 +98,42 @@ end @test successful_retcode(model) end +@testitem "leaky bucket model" begin + using SciMLBase: successful_retcode + import BasicModelInterface as BMI + + toml_path = normpath(@__DIR__, "../../generated_testmodels/leaky_bucket/ribasim.toml") + @test ispath(toml_path) + model = Ribasim.Model(toml_path) + @test model isa Ribasim.Model + + stor = model.integrator.u.storage + prec = model.integrator.p.basin.precipitation + evap = model.integrator.p.basin.potential_evaporation + drng = model.integrator.p.basin.drainage + infl = model.integrator.p.basin.infiltration + # The dynamic data has missings, but these are not set. + @test prec == [0.0] + @test evap == [0.0] + @test drng == [0.003] + @test infl == [0.0] + init_stor = 1000.0 + @test stor == [init_stor] + BMI.update_until(model, 1.5 * 86400) + @test prec == [0.0] + @test evap == [0.0] + @test drng == [0.003] + @test infl == [0.001] + stor ≈ Float32[init_stor + 86400 * (0.003 * 1.5 - 0.001 * 0.5)] + BMI.update_until(model, 2.5 * 86400) + @test prec == [0.00] + @test evap == [0.0] + @test drng == [0.001] + @test infl == [0.002] + stor ≈ Float32[init_stor + 86400 * (0.003 * 2.0 + 0.001 * 0.5 - 0.001 - 0.002 * 0.5)] + @test successful_retcode(Ribasim.solve!(model)) +end + @testitem "basic model" begin using Logging: Debug, with_logger using LoggingExtras diff --git a/docs/schema/BasinTime.schema.json b/docs/schema/BasinTime.schema.json index badcdf567..bfbcb03e2 100644 --- a/docs/schema/BasinTime.schema.json +++ b/docs/schema/BasinTime.schema.json @@ -14,24 +14,59 @@ "type": "string" }, "drainage": { - "format": "double", - "type": "number" + "format": "default", + "anyOf": [ + { + "type": "null" + }, + { + "type": "number" + } + ] }, "potential_evaporation": { - "format": "double", - "type": "number" + "format": "default", + "anyOf": [ + { + "type": "null" + }, + { + "type": "number" + } + ] }, "infiltration": { - "format": "double", - "type": "number" + "format": "default", + "anyOf": [ + { + "type": "null" + }, + { + "type": "number" + } + ] }, "precipitation": { - "format": "double", - "type": "number" + "format": "default", + "anyOf": [ + { + "type": "null" + }, + { + "type": "number" + } + ] }, "urban_runoff": { - "format": "double", - "type": "number" + "format": "default", + "anyOf": [ + { + "type": "null" + }, + { + "type": "number" + } + ] }, "remarks": { "description": "a hack for pandera", @@ -42,11 +77,6 @@ }, "required": [ "node_id", - "time", - "drainage", - "potential_evaporation", - "infiltration", - "precipitation", - "urban_runoff" + "time" ] } diff --git a/python/ribasim/ribasim/models.py b/python/ribasim/ribasim/models.py index 87605768b..f5c63246d 100644 --- a/python/ribasim/ribasim/models.py +++ b/python/ribasim/ribasim/models.py @@ -44,11 +44,11 @@ class BasinSubgrid(BaseModel): class BasinTime(BaseModel): node_id: int time: datetime - drainage: float - potential_evaporation: float - infiltration: float - precipitation: float - urban_runoff: float + drainage: float | None = None + potential_evaporation: float | None = None + infiltration: float | None = None + precipitation: float | None = None + urban_runoff: float | None = None remarks: str = Field("", description="a hack for pandera") diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index f8c1f7d35..a73777f05 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -21,7 +21,7 @@ outlet_model, tabulated_rating_curve_model, ) -from ribasim_testmodels.bucket import bucket_model +from ribasim_testmodels.bucket import bucket_model, leaky_bucket_model from ribasim_testmodels.discrete_control import ( flow_condition_model, level_boundary_condition_model, @@ -54,37 +54,38 @@ __all__ = [ "allocation_example_model", "backwater_model", - "basic_model", "basic_arrow_model", + "basic_model", "basic_transient_model", "bucket_model", - "pump_discrete_control_model", - "flow_condition_model", - "tabulated_rating_curve_model", - "trivial_model", - "linear_resistance_model", - "rating_curve_model", - "manning_resistance_model", - "pid_control_model", - "misc_nodes_model", - "tabulated_rating_curve_control_model", + "discrete_control_of_pid_control_model", "dutch_waterways_model", - "invalid_qh_model", "flow_boundary_time_model", - "pid_control_equation_model", - "invalid_fractional_flow_model", + "flow_condition_model", + "fractional_flow_subnetwork_model", "invalid_discrete_control_model", - "level_setpoint_with_minmax_model", "invalid_edge_types_model", - "discrete_control_of_pid_control_model", + "invalid_fractional_flow_model", + "invalid_qh_model", + "leaky_bucket_model", "level_boundary_condition_model", + "level_setpoint_with_minmax_model", + "linear_resistance_model", + "looped_subnetwork_model", + "manning_resistance_model", + "minimal_subnetwork_model", + "misc_nodes_model", "outlet_model", - "user_model", + "pid_control_equation_model", + "pid_control_model", + "pump_discrete_control_model", + "rating_curve_model", "subnetwork_model", - "minimal_subnetwork_model", - "fractional_flow_subnetwork_model", - "looped_subnetwork_model", + "tabulated_rating_curve_control_model", + "tabulated_rating_curve_model", + "trivial_model", "two_basin_model", + "user_model", ] # provide a mapping from model name to its constructor, so we can iterate over all models diff --git a/python/ribasim_testmodels/ribasim_testmodels/bucket.py b/python/ribasim_testmodels/ribasim_testmodels/bucket.py index 1dd5dce06..f16e06319 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/bucket.py +++ b/python/ribasim_testmodels/ribasim_testmodels/bucket.py @@ -76,3 +76,78 @@ def bucket_model() -> ribasim.Model: endtime="2021-01-01 00:00:00", ) return model + + +def leaky_bucket_model() -> ribasim.Model: + """Bucket model with dynamic forcing with missings.""" + + # Set up the nodes: + xy = np.array( + [ + (400.0, 200.0), # Basin + ] + ) + node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1]) + node_type = ["Basin"] + # Make sure the feature id starts at 1: explicitly give an index. + node = ribasim.Node( + df=gpd.GeoDataFrame( + data={"type": node_type}, + index=pd.Index(np.arange(len(xy)) + 1, name="fid"), + geometry=node_xy, + crs="EPSG:28992", + ) + ) + + # Setup the dummy edges: + from_id = np.array([], dtype=np.int64) + to_id = np.array([], dtype=np.int64) + lines = node.geometry_from_connectivity(from_id.tolist(), to_id.tolist()) + edge = ribasim.Edge( + df=gpd.GeoDataFrame( + data={ + "from_node_id": from_id, + "to_node_id": to_id, + "edge_type": len(from_id) * ["flow"], + }, + geometry=lines, + crs="EPSG:28992", + ) + ) + + # Setup the basins: + profile = pd.DataFrame( + data={ + "node_id": [1, 1], + "area": [1000.0, 1000.0], + "level": [0.0, 1.0], + } + ) + + state = pd.DataFrame( + data={ + "node_id": [1], + "level": [1.0], + } + ) + + time = pd.DataFrame( + data={ + "time": pd.date_range("2020-01-01", "2020-01-05"), + "node_id": 1, + "drainage": [0.003, np.nan, 0.001, 0.002, 0.0], + "potential_evaporation": np.nan, + "infiltration": [np.nan, 0.001, 0.002, 0.0, 0.0], + "precipitation": np.nan, + "urban_runoff": 0.0, + } + ) + basin = ribasim.Basin(profile=profile, time=time, state=state) + + model = ribasim.Model( + network=ribasim.Network(node=node, edge=edge), + basin=basin, + starttime="2020-01-01 00:00:00", + endtime="2020-01-05 00:00:00", + ) + return model From f721804eef3176a184e42147a7e1f54b1c8015ae Mon Sep 17 00:00:00 2001 From: Bart de Koning <74617371+SouthEndMusic@users.noreply.github.com> Date: Fri, 2 Feb 2024 13:27:06 +0100 Subject: [PATCH 3/4] Merge main allocation network feature branch (#1006) Fixes https://github.com/Deltares/Ribasim/issues/922. --------- Co-authored-by: Marnix Kraus Co-authored-by: Martijn Visser --- Manifest.toml | 2 +- core/src/allocation.jl | 393 +++++++--- core/src/bmi.jl | 43 +- core/src/create.jl | 32 +- core/src/solve.jl | 18 +- core/src/utils.jl | 15 + core/test/allocation_test.jl | 117 ++- docs/core/allocation.qmd | 6 +- docs/python/test-models.qmd | 2 +- python/ribasim/ribasim/geometry/node.py | 43 ++ python/ribasim/ribasim/model.py | 17 +- .../ribasim_testmodels/__init__.py | 2 + .../ribasim_testmodels/allocation.py | 691 +++++++++++++++++- 13 files changed, 1232 insertions(+), 149 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 4b4538c67..42cbf3109 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.0" manifest_format = "2.0" -project_hash = "7e18cea33bf679231c5c94c6b99a1b635f9a1b61" +project_hash = "742e16f5de5af91dfff30d91944981a5c7673d44" [[deps.ADTypes]] git-tree-sha1 = "41c37aa88889c171f1300ceac1313c06e891d245" diff --git a/core/src/allocation.jl b/core/src/allocation.jl index ac6a0745f..f8248bbc3 100644 --- a/core/src/allocation.jl +++ b/core/src/allocation.jl @@ -1,19 +1,45 @@ +"""Find the edges from the main network to a subnetwork.""" +function find_subnetwork_connections!(p::Parameters)::Nothing + (; allocation, graph, user) = p + n_priorities = length(user.demand[1]) + (; subnetwork_demands, subnetwork_allocateds) = allocation + for node_id in graph[].node_ids[1] + for outflow_id in outflow_ids(graph, node_id) + if graph[outflow_id].allocation_network_id != 1 + main_network_source_edges = + get_main_network_connections(p, graph[outflow_id].allocation_network_id) + edge = (node_id, outflow_id) + push!(main_network_source_edges, edge) + subnetwork_demands[edge] = zeros(n_priorities) + subnetwork_allocateds[edge] = zeros(n_priorities) + end + end + end + return nothing +end + """ Find all nodes in the subnetwork which will be used in the allocation network. Some nodes are skipped to optimize allocation optimization. """ function allocation_graph_used_nodes!(p::Parameters, allocation_network_id::Int)::Nothing - (; graph, basin, fractional_flow) = p + (; graph, basin, fractional_flow, allocation) = p + (; main_network_connections) = allocation node_ids = graph[].node_ids[allocation_network_id] used_nodes = Set{NodeID}() for node_id in node_ids + use_node = false has_fractional_flow_outneighbors = get_fractional_flow_connected_basins(node_id, basin, fractional_flow, graph)[3] node_type = graph[node_id].type if node_type in [:user, :basin, :terminal] - push!(used_nodes, node_id) + use_node = true elseif has_fractional_flow_outneighbors + use_node = true + end + + if use_node push!(used_nodes, node_id) end end @@ -21,13 +47,24 @@ function allocation_graph_used_nodes!(p::Parameters, allocation_network_id::Int) # Add nodes in the allocation graph for nodes connected to the source edges # One of these nodes can be outside the subnetwork, as long as the edge # connects to the subnetwork - for edge_metadata in graph[].edges_source[allocation_network_id] + edges_source = graph[].edges_source + for edge_metadata in get(edges_source, allocation_network_id, Set{EdgeMetadata}()) (; from_id, to_id) = edge_metadata push!(used_nodes, from_id) push!(used_nodes, to_id) end filter!(in(used_nodes), node_ids) + + # For the main network, include nodes that connect the main network to a subnetwork + # (also includes nodes not in the main network in the input) + if is_main_network(allocation_network_id) + for connections_subnetwork in main_network_connections + for connection in connections_subnetwork + union!(node_ids, connection) + end + end + end return nothing end @@ -241,6 +278,8 @@ function process_allocation_graph_edges!( return capacity end +const allocation_source_nodetypes = Set{Symbol}([:level_boundary, :flow_boundary]) + """ The source nodes must only have one allocation outneighbor and no allocation inneighbors. """ @@ -251,24 +290,23 @@ function valid_sources(p::Parameters, allocation_network_id::Int)::Bool errors = false - for (id_source, id_dst) in edge_ids + for edge in edge_ids + (id_source, id_dst) = edge if graph[id_source, id_dst].allocation_network_id_source == allocation_network_id - ids_allocation_in = [ - label for label in inneighbor_labels(graph, id_source) if - graph[label, id_source].allocation_flow - ] - if length(ids_allocation_in) !== 0 - errors = true - @error "Source edge ($id_source, $id_dst) is not an entry point of subnetwork $allocation_network_id" - end + from_source_node = graph[id_source].type in allocation_source_nodetypes - ids_allocation_out = [ - label for label in outneighbor_labels(graph, id_source) if - graph[id_source, label].allocation_flow - ] - if length(ids_allocation_out) !== 1 - errors = true - @error "Source edge ($id_source, $id_dst) is not the only allocation edge coming from $id_source" + if is_main_network(allocation_network_id) + if !from_source_node + errors = true + @error "The source node of source edge $edge in the main network must be one of $allocation_source_nodetypes." + end + else + from_main_network = is_main_network(graph[id_source].allocation_network_id) + + if !from_source_node && !from_main_network + errors = true + @error "The source node of source edge $edge for subnetwork $allocation_network_id is neither a source node nor is it coming from the main network." + end end end end @@ -297,6 +335,25 @@ function avoid_using_own_returnflow!(p::Parameters, allocation_network_id::Int): return nothing end +""" +Add the edges connecting the main network work to a subnetwork to both the main network +and subnetwork allocation graph. +""" +function add_subnetwork_connections!(p::Parameters, allocation_network_id::Int)::Nothing + (; graph, allocation) = p + (; main_network_connections) = allocation + edge_ids = graph[].edge_ids[allocation_network_id] + + if is_main_network(allocation_network_id) + for connections in main_network_connections + union!(edge_ids, connections) + end + else + union!(edge_ids, get_main_network_connections(p, allocation_network_id)) + end + return nothing +end + """ Build the graph used for the allocation problem. """ @@ -312,6 +369,7 @@ function allocation_graph( # Process the edges in the allocation graph process_allocation_graph_edges!(capacity, edges_composite, p, allocation_network_id) + add_subnetwork_connections!(p, allocation_network_id) if !valid_sources(p, allocation_network_id) error("Errors in sources in allocation graph.") @@ -351,10 +409,21 @@ function add_variables_absolute_value!( allocation_network_id::Int, config::Config, )::Nothing - (; graph) = p - node_ids = graph[].node_ids[allocation_network_id] - node_ids_user = [node_id for node_id in node_ids if graph[node_id].type == :user] + (; graph, allocation) = p + (; main_network_connections) = allocation if startswith(config.allocation.objective_type, "linear") + node_ids = graph[].node_ids[allocation_network_id] + node_ids_user = [node_id for node_id in node_ids if graph[node_id].type == :user] + + # For the main network, connections to subnetworks are treated as users + if is_main_network(allocation_network_id) + for connections_subnetwork in main_network_connections + for connection in connections_subnetwork + push!(node_ids_user, connection[2]) + end + end + end + problem[:F_abs] = JuMP.@variable(problem, F_abs[node_id = node_ids_user]) end return nothing @@ -375,11 +444,12 @@ function add_constraints_capacity!( allocation_network_id::Int, )::Nothing (; graph) = p + main_network_source_edges = get_main_network_connections(p, allocation_network_id) F = problem[:F] edge_ids = graph[].edge_ids[allocation_network_id] edge_ids_finite_capacity = Tuple{NodeID, NodeID}[] for edge in edge_ids - if !isinf(capacity[edge...]) + if !isinf(capacity[edge...]) && edge ∉ main_network_source_edges push!(edge_ids_finite_capacity, edge) end end @@ -461,13 +531,19 @@ function add_constraints_flow_conservation!( p::Parameters, allocation_network_id::Int, )::Nothing - (; graph) = p + (; graph, allocation) = p F = problem[:F] node_ids = graph[].node_ids[allocation_network_id] - node_ids_basin = [node_id for node_id in node_ids if graph[node_id].type == :basin] + node_ids_conservation = + [node_id for node_id in node_ids if graph[node_id].type == :basin] + main_network_source_edges = get_main_network_connections(p, allocation_network_id) + for edge in main_network_source_edges + push!(node_ids_conservation, edge[2]) + end + unique!(node_ids_conservation) problem[:flow_conservation] = JuMP.@constraint( problem, - [node_id = node_ids_basin], + [node_id = node_ids_conservation], sum([ F[(node_id, outneighbor_id)] for outneighbor_id in outflow_ids_allocation(graph, node_id) @@ -523,12 +599,23 @@ function add_constraints_absolute_value!( allocation_network_id::Int, config::Config, )::Nothing - (; graph) = p - node_ids = graph[].node_ids[allocation_network_id] + (; graph, allocation) = p + (; main_network_connections) = allocation objective_type = config.allocation.objective_type if startswith(objective_type, "linear") + node_ids = graph[].node_ids[allocation_network_id] node_ids_user = [node_id for node_id in node_ids if graph[node_id].type == :user] + + # For the main network, connections to subnetworks are treated as users + if is_main_network(allocation_network_id) + for connections_subnetwork in main_network_connections + for connection in connections_subnetwork + push!(node_ids_user, connection[2]) + end + end + end + node_ids_user_inflow = Dict( node_id_user => only(inflow_ids_allocation(graph, node_id_user)) for node_id_user in node_ids_user @@ -689,6 +776,60 @@ function AllocationModel( ) end +""" +Add a term to the expression of the objective function corresponding to +the demand of a user. +""" +function add_user_term!( + ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, + edge::Tuple{NodeID, NodeID}, + objective_type::Symbol, + demand::Float64, + model::AllocationModel, +)::Nothing + (; problem) = model + F = problem[:F] + F_edge = F[edge] + node_id_user = edge[2] + + if objective_type == :quadratic_absolute + # Objective function ∑ (F - d)^2 + JuMP.add_to_expression!(ex, 1, F_edge, F_edge) + JuMP.add_to_expression!(ex, -2 * demand, F_edge) + JuMP.add_to_expression!(ex, demand^2) + + elseif objective_type == :quadratic_relative + # Objective function ∑ (1 - F/d)^2 + if demand ≈ 0 + return nothing + end + JuMP.add_to_expression!(ex, 1.0 / demand^2, F_edge, F_edge) + JuMP.add_to_expression!(ex, -2.0 / demand, F_edge) + JuMP.add_to_expression!(ex, 1.0) + + elseif objective_type == :linear_absolute + # Objective function ∑ |F - d| + JuMP.set_normalized_rhs(problem[:abs_positive][node_id_user], -demand) + JuMP.set_normalized_rhs(problem[:abs_negative][node_id_user], demand) + + elseif objective_type == :linear_relative + # Objective function ∑ |1 - F/d| + JuMP.set_normalized_coefficient( + problem[:abs_positive][node_id_user], + F_edge, + iszero(demand) ? 0 : 1 / demand, + ) + JuMP.set_normalized_coefficient( + problem[:abs_negative][node_id_user], + F_edge, + iszero(demand) ? 0 : -1 / demand, + ) + else + error("Invalid allocation objective type $objective_type.") + end + return nothing +end + """ Set the objective for the given priority. For an objective with absolute values this also involves adjusting constraints. @@ -700,8 +841,9 @@ function set_objective_priority!( priority_idx::Int, )::Nothing (; objective_type, problem, allocation_network_id) = allocation_model - (; graph, user) = p + (; graph, user, allocation) = p (; demand, node_id) = user + (; main_network_connections, subnetwork_demands) = allocation edge_ids = graph[].edge_ids[allocation_network_id] F = problem[:F] @@ -713,6 +855,18 @@ function set_objective_priority!( demand_max = 0.0 + # Terms for subnetworks as users + if is_main_network(allocation_network_id) + for connections_subnetwork in main_network_connections + for connection in connections_subnetwork + d = subnetwork_demands[connection][priority_idx] + demand_max = max(demand_max, d) + add_user_term!(ex, connection, objective_type, d, allocation_model) + end + end + end + + # Terms for user nodes for edge_id in edge_ids node_id_user = edge_id[2] if graph[node_id_user].type != :user @@ -722,43 +876,7 @@ function set_objective_priority!( user_idx = findsorted(node_id, node_id_user) d = demand[user_idx][priority_idx](t) demand_max = max(demand_max, d) - F_edge = F[edge_id] - - if objective_type == :quadratic_absolute - # Objective function ∑ (F - d)^2 - JuMP.add_to_expression!(ex, 1, F_edge, F_edge) - JuMP.add_to_expression!(ex, -2 * d, F_edge) - JuMP.add_to_expression!(ex, d^2) - - elseif objective_type == :quadratic_relative - # Objective function ∑ (1 - F/d)^2S - if d ≈ 0 - continue - end - JuMP.add_to_expression!(ex, 1.0 / d^2, F_edge, F_edge) - JuMP.add_to_expression!(ex, -2.0 / d, F_edge) - JuMP.add_to_expression!(ex, 1.0) - - elseif objective_type == :linear_absolute - # Objective function ∑ |F - d| - JuMP.set_normalized_rhs(problem[:abs_positive][node_id_user], -d) - JuMP.set_normalized_rhs(problem[:abs_negative][node_id_user], d) - - elseif objective_type == :linear_relative - # Objective function ∑ |1 - F/d| - JuMP.set_normalized_coefficient( - problem[:abs_positive][node_id_user], - F_edge, - iszero(d) ? 0 : 1 / d, - ) - JuMP.set_normalized_coefficient( - problem[:abs_negative][node_id_user], - F_edge, - iszero(d) ? 0 : -1 / d, - ) - else - error("Invalid allocation objective type $objective_type.") - end + add_user_term!(ex, edge_id, objective_type, d, allocation_model) end # Add flow cost @@ -788,35 +906,66 @@ function assign_allocations!( allocation_model::AllocationModel, p::Parameters, t::Float64, - priority_idx::Int, + priority_idx::Int; + collect_demands::Bool = false, )::Nothing (; problem, allocation_network_id) = allocation_model - (; graph, user) = p + (; graph, user, allocation) = p + (; + subnetwork_demands, + subnetwork_allocateds, + allocation_network_ids, + main_network_connections, + ) = allocation (; record) = user edge_ids = graph[].edge_ids[allocation_network_id] + main_network_source_edges = get_main_network_connections(p, allocation_network_id) F = problem[:F] for edge_id in edge_ids + # If this edge is a source edge from the main network to a subnetwork, + # and demands are being collected, add its flow to the demand of this edge + if collect_demands && + graph[edge_id...].allocation_network_id_source == allocation_network_id && + edge_id ∈ main_network_source_edges + allocated = JuMP.value(F[edge_id]) + subnetwork_demands[edge_id][priority_idx] += allocated + end + user_node_id = edge_id[2] - if graph[user_node_id].type != :user - continue + + if graph[user_node_id].type == :user + allocated = JuMP.value(F[edge_id]) + user_idx = findsorted(user.node_id, user_node_id) + user.allocated[user_idx][priority_idx] = allocated + + # Save allocations to record + push!(record.time, t) + push!(record.allocation_network_id, allocation_model.allocation_network_id) + push!(record.user_node_id, Int(user_node_id)) + push!(record.priority, user.priorities[priority_idx]) + push!(record.demand, user.demand[user_idx][priority_idx](t)) + push!(record.allocated, allocated) + # TODO: This is now the last abstraction before the allocation update, + # should be the average abstraction since the last allocation solve + push!( + record.abstracted, + get_flow(graph, inflow_id(graph, user_node_id), user_node_id, 0), + ) + end + end + + # Write the flows to the subnetworks as allocated flows + # in the allocation object + if is_main_network(allocation_network_id) + for (allocation_network_id, main_network_source_edges) in + zip(allocation_network_ids, main_network_connections) + if is_main_network(allocation_network_id) + continue + end + for edge_id in main_network_source_edges + subnetwork_allocateds[edge_id][priority_idx] = JuMP.value(F[edge_id]) + end end - user_idx = findsorted(user.node_id, user_node_id) - allocated = JuMP.value(F[edge_id]) - user.allocated[user_idx][priority_idx] = allocated - - # Save allocations to record - push!(record.time, t) - push!(record.allocation_network_id, allocation_model.allocation_network_id) - push!(record.user_node_id, Int(user_node_id)) - push!(record.priority, user.priorities[priority_idx]) - push!(record.demand, user.demand[user_idx][priority_idx](t)) - push!(record.allocated, allocated) - # TODO: This is now the last abstraction before the allocation update, - # should be the average abstraction since the last allocation solve - push!( - record.abstracted, - get_flow(graph, inflow_id(graph, user_node_id), user_node_id, 0), - ) end return nothing end @@ -824,27 +973,43 @@ end """ Adjust the source flows. """ -function adjust_source_flows!( +function adjust_source_capacities!( allocation_model::AllocationModel, p::Parameters, - priority_idx::Int, + priority_idx::Int; + collect_demands::Bool = false, )::Nothing (; problem) = allocation_model - (; graph) = p + (; graph, allocation) = p (; allocation_network_id) = allocation_model + (; subnetwork_allocateds) = allocation edge_ids = graph[].edge_ids[allocation_network_id] source_constraints = problem[:source] F = problem[:F] - # It is assumed that the allocation procedure does not have to be differentiated. + main_network_source_edges = get_main_network_connections(p, allocation_network_id) + for edge_id in edge_ids - # If it is a source edge. if graph[edge_id...].allocation_network_id_source == allocation_network_id + # If it is a source edge for this allocation problem if priority_idx == 1 - # Reset the source to the current flow. + # If the optimization was just started, i.e. sources have to be reset + if edge_id in main_network_source_edges + if collect_demands + # Set the source capacity to effectively unlimited if subnetwork demands are being collected + source_capacity = Inf + else + # Set the source capacity to the value allocated to the subnetwork over this edge + source_capacity = subnetwork_allocateds[edge_id][priority_idx] + end + else + # Reset the source to the current flow from the physical layer. + source_capacity = get_flow(graph, edge_id..., 0) + end JuMP.set_normalized_rhs( source_constraints[edge_id], - get_flow(graph, edge_id..., 0), + # It is assumed that the allocation procedure does not have to be differentiated. + source_capacity, ) else # Subtract the allocated flow from the source. @@ -876,11 +1041,15 @@ function adjust_edge_capacities!( constraints_capacity = problem[:capacity] F = problem[:F] + main_network_source_edges = get_main_network_connections(p, allocation_network_id) + for edge_id in edge_ids c = capacity[edge_id...] - # Edges with infinite capacity have no capacity constraints - if isinf(c) + # These edges have no capacity constraints: + # - With infinite capacity + # - Being a source from the main network to a subnetwork + if isinf(c) || edge_id ∈ main_network_source_edges continue end @@ -902,18 +1071,34 @@ 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 - (; user) = p - (; problem) = allocation_model +function allocate!( + p::Parameters, + allocation_model::AllocationModel, + t::Float64; + collect_demands::Bool = false, +)::Nothing + (; user, allocation) = p + (; problem, allocation_network_id) = allocation_model (; priorities) = user + (; subnetwork_demands) = allocation # TODO: Compute basin flow from vertical fluxes and basin volume. # Set as basin demand if the net flow is negative, set as source # in the flow_conservation constraints if the net flow is positive. # Solve this as a separate problem before the priorities below + main_network_source_edges = get_main_network_connections(p, allocation_network_id) + + if collect_demands + for main_network_connection in keys(subnetwork_demands) + if main_network_connection in main_network_source_edges + subnetwork_demands[main_network_connection] .= 0.0 + end + end + end + for priority_idx in eachindex(priorities) - adjust_source_flows!(allocation_model, p, priority_idx) + adjust_source_capacities!(allocation_model, p, priority_idx; collect_demands) # Subtract the flows used by the allocation of the previous priority from the capacities of the edges # or set edge capacities if priority_idx = 1 @@ -930,10 +1115,14 @@ function allocate!(p::Parameters, allocation_model::AllocationModel, t::Float64) JuMP.optimize!(problem) @debug JuMP.solution_summary(problem) if JuMP.termination_status(problem) !== JuMP.OPTIMAL - error("Allocation coudn't find optimal solution.") + (; allocation_network_id) = allocation_model + priority = priorities[priority_index] + error( + "Allocation of subnetwork $allocation_network_id, priority $priority coudn't find optimal solution.", + ) end # Assign the allocations to the users for this priority - assign_allocations!(allocation_model, p, t, priority_idx) + assign_allocations!(allocation_model, p, t, priority_idx; collect_demands) end end diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 24c3b763a..c014dad0d 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -476,16 +476,30 @@ function discrete_control_affect!( return control_state_change end -function get_allocation_model( +function get_allocation_model(p::Parameters, allocation_network_id::Int)::AllocationModel + (; allocation) = p + (; allocation_network_ids, allocation_models) = allocation + idx = findsorted(allocation_network_ids, allocation_network_id) + if isnothing(idx) + error("Invalid allocation network ID $allocation_network_id.") + else + return allocation_models[idx] + end +end + +function get_main_network_connections( p::Parameters, allocation_network_id::Int, -)::Union{AllocationModel, Nothing} - for allocation_model in p.allocation_models - if allocation_model.allocation_network_id == allocation_network_id - return allocation_model - end +)::Vector{Tuple{NodeID, NodeID}} + (; allocation) = p + (; allocation_network_ids, main_network_connections) = allocation + idx = findsorted(allocation_network_ids, allocation_network_id) + if isnothing(idx) + error("Invalid allocation network ID $allocation_network_id.") + else + return main_network_connections[idx] end - return nothing + return end """ @@ -589,7 +603,20 @@ end "Solve the allocation problem for all users and assign allocated abstractions to user nodes." function update_allocation!(integrator)::Nothing (; p, t) = integrator - for allocation_model in integrator.p.allocation_models + (; allocation) = p + (; allocation_models) = allocation + + # If a main network is present, collect demands of subnetworks + if has_main_network(allocation) + for allocation_model in Iterators.drop(allocation_models, 1) + allocate!(p, allocation_model, t; collect_demands = true) + end + end + + # Solve the allocation problems + # If a main network is present this is solved first, + # which provides allocation to the subnetworks + for allocation_model in allocation_models allocate!(p, allocation_model, t) end end diff --git a/core/src/create.jl b/core/src/create.jl index 331a944f3..7c5494117 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -196,16 +196,26 @@ end const nonconservative_nodetypes = Set{String}(["Basin", "LevelBoundary", "FlowBoundary", "Terminal", "User"]) -function generate_allocation_models!(p::Parameters, config::Config)::Nothing - (; graph, allocation_models) = p +function initialize_allocation!(p::Parameters, config::Config)::Nothing + (; graph, allocation) = p + (; allocation_network_ids, allocation_models, main_network_connections) = allocation + allocation_network_ids_ = sort(collect(keys(graph[].node_ids))) errors = non_positive_allocation_network_id(graph) - if errors error("Allocation network initialization failed.") end - for allocation_network_id in keys(graph[].node_ids) + for allocation_network_id in allocation_network_ids_ + push!(allocation_network_ids, allocation_network_id) + push!(main_network_connections, Tuple{NodeID, NodeID}[]) + end + + if first(allocation_network_ids_) == 1 + find_subnetwork_connections!(p) + end + + for allocation_network_id in allocation_network_ids_ push!( allocation_models, AllocationModel(config, allocation_network_id, p, config.allocation.timestep), @@ -634,7 +644,7 @@ function User(db::DB, config::Config)::User error("Problems encountered when parsing User static and time node IDs.") end - # The highest priority number given, which corresponds to the least important demands + # All provided priorities priorities = sort(unique(union(static.priority, time.priority))) active = BitVector() @@ -798,7 +808,13 @@ 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_models = Vector{AllocationModel}() + allocation = Allocation( + Int[], + AllocationModel[], + Vector{Tuple{NodeID, NodeID}}[], + Dict{Tuple{NodeID, NodeID}, Float64}(), + Dict{Tuple{NodeID, NodeID}, Float64}(), + ) if !valid_edges(graph) error("Invalid edge(s) found.") @@ -835,7 +851,7 @@ function Parameters(db::DB, config::Config)::Parameters p = Parameters( config.starttime, graph, - allocation_models, + allocation, basin, linear_resistance, manning_resistance, @@ -859,7 +875,7 @@ function Parameters(db::DB, config::Config)::Parameters # Allocation data structures if config.allocation.use_allocation - generate_allocation_models!(p, config) + initialize_allocation!(p, config) end return p end diff --git a/core/src/solve.jl b/core/src/solve.jl index 2087f2f35..d766a86d5 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -21,6 +21,22 @@ struct AllocationModel Δt_allocation::Float64 end +""" +Object for all information about allocation +allocation_network_ids: The unique sorted allocation network IDs +allocation models: The allocation models for the main network and subnetworks corresponding to + allocation_network_ids +main_network_connections: (from_id, to_id) from the main network to the subnetwork per subnetwork +subnetwork_demands: The demand of an edge from the main network to a subnetwork +""" +struct Allocation + allocation_network_ids::Vector{Int} + allocation_models::Vector{AllocationModel} + main_network_connections::Vector{Vector{Tuple{NodeID, NodeID}}} + subnetwork_demands::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} + subnetwork_allocateds::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} +end + @enumx EdgeType flow control none """ @@ -471,7 +487,7 @@ struct Parameters{T, C1, C2} MetaGraphsNext.var"#11#13", Float64, } - allocation_models::Vector{AllocationModel} + allocation::Allocation basin::Basin{T, C1} linear_resistance::LinearResistance manning_resistance::ManningResistance diff --git a/core/src/utils.jl b/core/src/utils.jl index 27632f41b..f2f743da4 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -42,12 +42,19 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra db, "SELECT fid, from_node_id, to_node_id, edge_type, allocation_network_id FROM Edge ORDER BY fid", ) + # Node IDs per subnetwork node_ids = Dict{Int, Set{NodeID}}() + # Allocation edges per subnetwork edge_ids = Dict{Int, Set{Tuple{NodeID, NodeID}}}() + # Source edges per subnetwork edges_source = Dict{Int, Set{EdgeMetadata}}() + # The number of flow edges flow_counter = 0 + # Dictionary from flow edge to index in flow vector flow_dict = Dict{Tuple{NodeID, NodeID}, Int}() + # The number of nodes with vertical flow (interaction with outside of model) flow_vertical_counter = 0 + # Dictionary from node ID to index in vertical flow vector flow_vertical_dict = Dict{NodeID, Int}() graph = MetaGraph( DiGraph(); @@ -1260,3 +1267,11 @@ function allocation_path_exists_in_graph( end return false end + +function has_main_network(allocation::Allocation)::Bool + return first(allocation.allocation_network_ids) == 1 +end + +function is_main_network(allocation_network_id::Int)::Bool + return allocation_network_id == 1 +end diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index ff48730e6..5e1eaf0de 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -1,5 +1,4 @@ @testitem "Allocation solve" begin - using PreallocationTools: get_tmp using Ribasim: NodeID import SQLite import JuMP @@ -15,7 +14,7 @@ graph = p.graph Ribasim.set_flow!(graph, NodeID(1), NodeID(2), 4.5) # Source flow - allocation_model = p.allocation_models[1] + allocation_model = p.allocation.allocation_models[1] Ribasim.allocate!(p, allocation_model, 0.0) F = allocation_model.problem[:F] @@ -45,7 +44,7 @@ end config = Ribasim.Config(toml_path; allocation_objective_type = "quadratic_absolute") model = Ribasim.run(config) @test successful_retcode(model) - problem = model.integrator.p.allocation_models[1].problem + problem = model.integrator.p.allocation.allocation_models[1].problem objective = JuMP.objective_function(problem) @test objective isa JuMP.QuadExpr # Quadratic expression F = problem[:F] @@ -61,7 +60,7 @@ end config = Ribasim.Config(toml_path; allocation_objective_type = "quadratic_relative") model = Ribasim.run(config) @test successful_retcode(model) - problem = model.integrator.p.allocation_models[1].problem + problem = model.integrator.p.allocation.allocation_models[1].problem objective = JuMP.objective_function(problem) @test objective isa JuMP.QuadExpr # Quadratic expression @test objective.aff.constant == 2.0 @@ -78,7 +77,7 @@ end config = Ribasim.Config(toml_path; allocation_objective_type = "linear_absolute") model = Ribasim.run(config) @test successful_retcode(model) - problem = model.integrator.p.allocation_models[1].problem + problem = model.integrator.p.allocation.allocation_models[1].problem objective = JuMP.objective_function(problem) @test objective isa JuMP.AffExpr # Affine expression @test :F_abs in keys(problem.obj_dict) @@ -95,7 +94,7 @@ end config = Ribasim.Config(toml_path; allocation_objective_type = "linear_relative") model = Ribasim.run(config) @test successful_retcode(model) - problem = model.integrator.p.allocation_models[1].problem + problem = model.integrator.p.allocation.allocation_models[1].problem objective = JuMP.objective_function(problem) @test objective isa JuMP.AffExpr # Affine expression @test :F_abs in keys(problem.obj_dict) @@ -121,7 +120,7 @@ end "../../generated_testmodels/fractional_flow_subnetwork/ribasim.toml", ) model = Ribasim.BMI.initialize(Ribasim.Model, toml_path) - problem = model.integrator.p.allocation_models[1].problem + problem = model.integrator.p.allocation.allocation_models[1].problem F = problem[:F] @test JuMP.normalized_coefficient( problem[:fractional_flow][(NodeID(3), NodeID(5))], @@ -155,7 +154,7 @@ end @test record_control.control_state == ["A", "B"] fractional_flow_constraints = - model.integrator.p.allocation_models[1].problem[:fractional_flow] + model.integrator.p.allocation.allocation_models[1].problem[:fractional_flow] @test JuMP.normalized_coefficient( problem[:fractional_flow][(NodeID(3), NodeID(5))], F[(NodeID(2), NodeID(3))], @@ -165,3 +164,105 @@ end F[(NodeID(2), NodeID(3))], ) ≈ -0.25 end + +@testitem "main allocation network initialization" begin + using SQLite + using Ribasim: NodeID + + toml_path = normpath( + @__DIR__, + "../../generated_testmodels/main_network_with_subnetworks/ribasim.toml", + ) + @test ispath(toml_path) + cfg = Ribasim.Config(toml_path) + db_path = Ribasim.input_path(cfg, cfg.database) + db = SQLite.DB(db_path) + p = Ribasim.Parameters(db, cfg) + close(db) + (; allocation, graph) = p + (; main_network_connections, allocation_network_ids) = allocation + @test Ribasim.has_main_network(allocation) + @test Ribasim.is_main_network(first(allocation_network_ids)) + + # Connections from main network to subnetworks + @test isempty(main_network_connections[1]) + @test only(main_network_connections[2]) == (NodeID(2), NodeID(11)) + @test only(main_network_connections[3]) == (NodeID(6), NodeID(24)) + @test only(main_network_connections[4]) == (NodeID(10), NodeID(38)) + + # main-sub connections are part of main network allocation graph + allocation_edges_main_network = graph[].edge_ids[1] + @test Tuple{NodeID, NodeID}[(2, 11), (6, 24), (10, 38)] ⊆ allocation_edges_main_network + + # Subnetworks interpreted as users require variables and constraints to + # support absolute value expressions in the objective function + allocation_model_main_network = Ribasim.get_allocation_model(p, 1) + problem = allocation_model_main_network.problem + @test problem[:F_abs].axes[1] == NodeID[11, 24, 38] + @test problem[:abs_positive].axes[1] == NodeID[11, 24, 38] + @test problem[:abs_negative].axes[1] == NodeID[11, 24, 38] + + # In each subnetwork, the connection from the main network to the subnetwork is + # interpreted as a source + @test Ribasim.get_allocation_model(p, 3).problem[:source].axes[1] == + Tuple{NodeID, NodeID}[(2, 11)] + @test Ribasim.get_allocation_model(p, 5).problem[:source].axes[1] == + Tuple{NodeID, NodeID}[(6, 24)] + @test Ribasim.get_allocation_model(p, 7).problem[:source].axes[1] == + Tuple{NodeID, NodeID}[(10, 38)] +end + +@testitem "allocation with main network optimization problem" begin + using SQLite + using Ribasim: NodeID + using JuMP + + toml_path = normpath( + @__DIR__, + "../../generated_testmodels/main_network_with_subnetworks/ribasim.toml", + ) + @test ispath(toml_path) + cfg = Ribasim.Config(toml_path) + db_path = Ribasim.input_path(cfg, cfg.database) + db = SQLite.DB(db_path) + p = Ribasim.Parameters(db, cfg) + close(db) + + (; allocation, user, graph) = p + (; allocation_models, subnetwork_demands, subnetwork_allocateds) = allocation + t = 0.0 + + # Collecting demands + for allocation_model in allocation_models[2:end] + Ribasim.allocate!(p, allocation_model, t; collect_demands = true) + end + + @test subnetwork_demands[(NodeID(2), NodeID(11))] ≈ [4.0, 4.0, 0.0] + @test subnetwork_demands[(NodeID(6), NodeID(24))] ≈ [0.001333333333, 0.0, 0.0] + @test subnetwork_demands[(NodeID(10), NodeID(38))] ≈ [0.001, 0.002, 0.002] + + # Solving for the main network, + # containing subnetworks as users + allocation_model = allocation_models[1] + (; problem) = allocation_model + Ribasim.allocate!(p, allocation_model, t) + + # Main network objective function + objective = JuMP.objective_function(problem) + objective_variables = keys(objective.terms) + F_abs = problem[:F_abs] + @test F_abs[NodeID(11)] ∈ objective_variables + @test F_abs[NodeID(24)] ∈ objective_variables + @test F_abs[NodeID(38)] ∈ objective_variables + + # Running full allocation algorithm + Ribasim.set_flow!(graph, NodeID(1), NodeID(2), 4.5) + Ribasim.update_allocation!((; p, t)) + + @test subnetwork_allocateds[NodeID(2), NodeID(11)] ≈ [4.0, 0.49766666, 0.0] + @test subnetwork_allocateds[NodeID(6), NodeID(24)] ≈ [0.00133333333, 0.0, 0.0] + @test subnetwork_allocateds[NodeID(10), NodeID(38)] ≈ [0.001, 0.0, 0.0] + + @test user.allocated[2] ≈ [4.0, 0.0, 0.0] + @test user.allocated[7] ≈ [0.001, 0.0, 0.0] +end diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index 14b16167e..e09d11765 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -238,15 +238,15 @@ using SQLite toml_path = normpath(@__DIR__, "../../generated_testmodels/allocation_example/ribasim.toml") p = Ribasim.Model(toml_path).integrator.p -allocation_model = p.allocation_models[1] +allocation_model = p.allocation.allocation_models[1] t = 0.0 priority_idx = 1 Ribasim.set_flow!(p.graph, Ribasim.NodeID(1), Ribasim.NodeID(2), 1.0) -Ribasim.adjust_source_flows!(allocation_model, p, priority_idx) +Ribasim.adjust_source_capacities!(allocation_model, p, priority_idx) Ribasim.adjust_edge_capacities!(allocation_model, p, priority_idx) Ribasim.set_objective_priority!(allocation_model, p, t, priority_idx) -println(p.allocation_models[1].problem) +println(p.allocation.allocation_models[1].problem) ``` diff --git a/docs/python/test-models.qmd b/docs/python/test-models.qmd index 6a60503cf..81891cd1d 100644 --- a/docs/python/test-models.qmd +++ b/docs/python/test-models.qmd @@ -2,7 +2,7 @@ title: "Test models" --- -Ribasim developers use the following models in its testbench and in order to test new features. +Ribasim developers use the following models in their testbench and in order to test new features. ```{python} # | code-fold: true diff --git a/python/ribasim/ribasim/geometry/node.py b/python/ribasim/ribasim/geometry/node.py index a4f9f552e..f10a2bc22 100644 --- a/python/ribasim/ribasim/geometry/node.py +++ b/python/ribasim/ribasim/geometry/node.py @@ -1,11 +1,13 @@ from collections.abc import Sequence from typing import Any +import geopandas as gpd import matplotlib.pyplot as plt import numpy as np import pandas as pd import pandera as pa import shapely +from matplotlib.patches import Patch from numpy.typing import NDArray from pandera.typing import Series from pandera.typing.geopandas import GeoSeries @@ -130,6 +132,47 @@ def connectivity_from_geometry( to_id = node_index[edge_node_id[:, 1]].to_numpy() return from_id, to_id + def plot_allocation_networks(self, ax=None, zorder=None) -> Any: + if ax is None: + _, ax = plt.subplots() + ax.axis("off") + + COLOR_SUBNETWORK = "black" + COLOR_MAIN_NETWORK = "blue" + ALPHA = 0.25 + + contains_main_network = False + contains_subnetworks = False + + for allocation_subnetwork_id, df_subnetwork in self.df.groupby( + "allocation_network_id" + ): + if allocation_subnetwork_id is None: + continue + elif allocation_subnetwork_id == 1: + contains_main_network = True + color = COLOR_MAIN_NETWORK + else: + contains_subnetworks = True + color = COLOR_SUBNETWORK + + hull = gpd.GeoDataFrame( + geometry=[df_subnetwork.geometry.unary_union.convex_hull] + ) + hull.plot(ax=ax, color=color, alpha=ALPHA, zorder=zorder) + + handles = [] + labels = [] + + if contains_main_network: + handles.append(Patch(facecolor=COLOR_MAIN_NETWORK, alpha=ALPHA)) + labels.append("Main network") + if contains_subnetworks: + handles.append(Patch(facecolor=COLOR_SUBNETWORK, alpha=ALPHA)) + labels.append("Subnetwork") + + return handles, labels + def plot(self, ax=None, zorder=None) -> Any: """ Plot the nodes. Each node type is given a separate marker. diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py index 5278b8efd..13b4509df 100644 --- a/python/ribasim/ribasim/model.py +++ b/python/ribasim/ribasim/model.py @@ -427,9 +427,9 @@ def plot_control_listen(self, ax): label="Listen edge" if i == 0 else None, ) - def plot(self, ax=None) -> Any: + def plot(self, ax=None, indicate_subnetworks: bool = True) -> Any: """ - Plot the nodes and edges of the model. + Plot the nodes, edges and allocation networks of the model. Parameters ---------- @@ -443,11 +443,22 @@ def plot(self, ax=None) -> Any: if ax is None: _, ax = plt.subplots() ax.axis("off") + self.network.edge.plot(ax=ax, zorder=2) self.plot_control_listen(ax) self.network.node.plot(ax=ax, zorder=3) - ax.legend(loc="lower left", bbox_to_anchor=(1, 0.5)) + handles, labels = ax.get_legend_handles_labels() + + if indicate_subnetworks: + ( + handles_subnetworks, + labels_subnetworks, + ) = self.network.node.plot_allocation_networks(ax=ax, zorder=1) + handles += handles_subnetworks + labels += labels_subnetworks + + ax.legend(handles, labels, loc="lower left", bbox_to_anchor=(1, 0.5)) return ax diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index a73777f05..4616d2746 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -9,6 +9,7 @@ allocation_example_model, fractional_flow_subnetwork_model, looped_subnetwork_model, + main_network_with_subnetworks_model, minimal_subnetwork_model, subnetwork_model, user_model, @@ -86,6 +87,7 @@ "trivial_model", "two_basin_model", "user_model", + "main_network_with_subnetworks_model", ] # provide a mapping from model name to its constructor, so we can iterate over all models diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 3be17b397..0e556823e 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -123,7 +123,9 @@ def user_model(): def subnetwork_model(): - """Create a user testmodel representing a subnetwork.""" + """Create a user testmodel representing a subnetwork. + This model is merged into main_network_with_subnetworks_model. + """ # Setup the nodes: xy = np.array( @@ -164,7 +166,7 @@ def subnetwork_model(): # Make sure the feature id starts at 1: explicitly give an index. node = ribasim.Node( df=gpd.GeoDataFrame( - data={"type": node_type, "allocation_network_id": 1}, + data={"type": node_type, "allocation_network_id": 2}, index=pd.Index(np.arange(len(xy)) + 1, name="fid"), geometry=node_xy, crs="EPSG:28992", @@ -177,7 +179,7 @@ def subnetwork_model(): ) to_id = np.array([2, 3, 4, 10, 5, 6, 7, 8, 11, 12, 13, 9, 2, 6, 8], dtype=np.int64) allocation_network_id = len(from_id) * [None] - allocation_network_id[0] = 1 + allocation_network_id[0] = 2 lines = node.geometry_from_connectivity(from_id, to_id) edge = ribasim.Edge( df=gpd.GeoDataFrame( @@ -283,7 +285,9 @@ def subnetwork_model(): def looped_subnetwork_model(): - """Create a user testmodel representing a subnetwork containing a loop in the topology.""" + """Create a user testmodel representing a subnetwork containing a loop in the topology. + This model is merged into main_network_with_subnetworks_model. + """ # Setup the nodes: xy = np.array( [ @@ -345,7 +349,7 @@ def looped_subnetwork_model(): # Make sure the feature id starts at 1: explicitly give an index. node = ribasim.Node( df=gpd.GeoDataFrame( - data={"type": node_type, "allocation_network_id": 1}, + data={"type": node_type, "allocation_network_id": 2}, index=pd.Index(np.arange(len(xy)) + 1, name="fid"), geometry=node_xy, crs="EPSG:28992", @@ -423,7 +427,7 @@ def looped_subnetwork_model(): ) lines = node.geometry_from_connectivity(from_id, to_id) allocation_network_id = len(from_id) * [None] - allocation_network_id[0] = 1 + allocation_network_id[0] = 2 edge = ribasim.Edge( df=gpd.GeoDataFrame( data={ @@ -556,7 +560,7 @@ def minimal_subnetwork_model(): # Make sure the feature id starts at 1: explicitly give an index. node = ribasim.Node( df=gpd.GeoDataFrame( - data={"type": node_type, "allocation_network_id": 1}, + data={"type": node_type, "allocation_network_id": 2}, index=pd.Index(np.arange(len(xy)) + 1, name="fid"), geometry=node_xy, crs="EPSG:28992", @@ -573,7 +577,7 @@ def minimal_subnetwork_model(): dtype=np.int64, ) allocation_network_id = len(from_id) * [None] - allocation_network_id[0] = 1 + allocation_network_id[0] = 2 lines = node.geometry_from_connectivity(from_id, to_id) edge = ribasim.Edge( df=gpd.GeoDataFrame( @@ -677,7 +681,9 @@ def minimal_subnetwork_model(): def fractional_flow_subnetwork_model(): - """Create a small subnetwork that contains fractional flow nodes.""" + """Create a small subnetwork that contains fractional flow nodes. + This model is merged into main_network_with_subnetworks_model. + """ xy = np.array( [ @@ -711,7 +717,7 @@ def fractional_flow_subnetwork_model(): # Make sure the feature id starts at 1: explicitly give an index. node = ribasim.Node( df=gpd.GeoDataFrame( - data={"type": node_type, "allocation_network_id": 1}, + data={"type": node_type, "allocation_network_id": 2}, index=pd.Index(np.arange(len(xy)) + 1, name="fid"), geometry=node_xy, crs="EPSG:28992", @@ -728,7 +734,7 @@ def fractional_flow_subnetwork_model(): dtype=np.int64, ) allocation_network_id = len(from_id) * [None] - allocation_network_id[0] = 1 + allocation_network_id[0] = 2 lines = node.geometry_from_connectivity(from_id, to_id) edge = ribasim.Edge( df=gpd.GeoDataFrame( @@ -903,10 +909,10 @@ def allocation_example_model(): "User", ] - # All nodes belong to allocation network id 1 + # All nodes belong to allocation network id 2 node = ribasim.Node( df=gpd.GeoDataFrame( - data={"type": node_type, "allocation_network_id": 1}, + data={"type": node_type, "allocation_network_id": 2}, index=pd.Index(np.arange(len(xy)) + 1, name="fid"), geometry=node_xy, crs="EPSG:28992", @@ -924,7 +930,7 @@ def allocation_example_model(): # Denote the first edge, 1 => 2, as a source edge for # allocation network 1 allocation_network_id = len(from_id) * [None] - allocation_network_id[0] = 1 + allocation_network_id[0] = 2 lines = node.geometry_from_connectivity(from_id, to_id) edge = ribasim.Edge( df=gpd.GeoDataFrame( @@ -1072,3 +1078,660 @@ def allocation_example_model(): ) return model + + +def main_network_with_subnetworks_model(): + # Set ip the nodes: + xy = np.array( + [ + (0.0, -1.0), + (3.0, 1.0), + (6.0, -1.0), + (9.0, 1.0), + (12.0, -1.0), + (15.0, 1.0), + (18.0, -1.0), + (21.0, 1.0), + (24.0, -1.0), + (27.0, 1.0), + (3.0, 4.0), + (2.0, 4.0), + (1.0, 4.0), + (0.0, 4.0), + (2.0, 5.0), + (2.0, 6.0), + (1.0, 6.0), + (0.0, 6.0), + (2.0, 8.0), + (2.0, 3.0), + (3.0, 6.0), + (0.0, 7.0), + (2.0, 7.0), + (14.0, 3.0), + (14.0, 4.0), + (14.0, 5.0), + (13.0, 6.0), + (12.0, 7.0), + (11.0, 8.0), + (15.0, 6.0), + (16.0, 7.0), + (17.0, 8.0), + (13.0, 5.0), + (26.0, 3.0), + (26.0, 4.0), + (25.0, 4.0), + (24.0, 4.0), + (28.0, 4.0), + (26.0, 5.0), + (24.0, 6.0), + (25.0, 6.0), + (26.0, 6.0), + (27.0, 6.0), + (28.0, 6.0), + (24.0, 7.0), + (26.0, 7.0), + (28.0, 7.0), + (26.0, 8.0), + (27.0, 8.0), + (28.0, 8.0), + (25.0, 9.0), + (26.0, 9.0), + (28.0, 9.0), + (26.0, 10.0), + (26.0, 11.0), + (26.0, 12.0), + (29.0, 6.0), + ] + ) + node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1]) + + node_type = [ + "FlowBoundary", + "Basin", + "LinearResistance", + "Basin", + "LinearResistance", + "Basin", + "LinearResistance", + "Basin", + "LinearResistance", + "Basin", + "Pump", + "Basin", + "Outlet", + "Terminal", + "Pump", + "Basin", + "Outlet", + "Basin", + "Terminal", + "User", + "User", + "User", + "Outlet", + "Pump", + "Basin", + "TabulatedRatingCurve", + "FractionalFlow", + "Basin", + "User", + "FractionalFlow", + "Basin", + "User", + "DiscreteControl", + "User", + "Basin", + "Outlet", + "Terminal", + "Pump", + "Pump", + "Basin", + "Outlet", + "Basin", + "Outlet", + "Basin", + "User", + "TabulatedRatingCurve", + "TabulatedRatingCurve", + "Basin", + "Pump", + "Basin", + "User", + "TabulatedRatingCurve", + "User", + "Basin", + "Outlet", + "Terminal", + "User", + ] + + allocation_network_id = np.ones(57, dtype=int) + allocation_network_id[10:23] = 3 + allocation_network_id[23:33] = 5 + allocation_network_id[33:] = 7 + + # Make sure the feature id starts at 1: explicitly give an index. + node = ribasim.Node( + df=gpd.GeoDataFrame( + data={ + "type": node_type, + "allocation_network_id": allocation_network_id, + }, + index=pd.Index(np.arange(len(xy)) + 1, name="fid"), + geometry=node_xy, + crs="EPSG:28992", + ) + ) + + # Setup the edges: + from_id = np.array( + [ + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 11, + 12, + 13, + 12, + 12, + 15, + 16, + 17, + 16, + 18, + 16, + 23, + 20, + 21, + 22, + 24, + 25, + 26, + 27, + 28, + 29, + 26, + 30, + 31, + 32, + 33, + 33, + 38, + 35, + 36, + 35, + 35, + 39, + 42, + 41, + 40, + 42, + 46, + 48, + 49, + 50, + 48, + 52, + 48, + 51, + 54, + 55, + 42, + 43, + 44, + 47, + 34, + 45, + 53, + 44, + 57, + 2, + 6, + 10, + ], + dtype=np.int64, + ) + to_id = np.array( + [ + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 12, + 13, + 14, + 20, + 15, + 16, + 17, + 18, + 21, + 22, + 23, + 19, + 12, + 16, + 18, + 25, + 26, + 27, + 28, + 29, + 28, + 30, + 31, + 32, + 31, + 27, + 30, + 35, + 36, + 37, + 34, + 39, + 42, + 41, + 40, + 45, + 46, + 48, + 49, + 50, + 53, + 52, + 54, + 51, + 54, + 55, + 56, + 43, + 44, + 47, + 50, + 35, + 40, + 50, + 57, + 44, + 11, + 24, + 38, + ], + dtype=np.int64, + ) + + edge_type = 68 * ["flow"] + edge_type[34] = "control" + edge_type[35] = "control" + allocation_network_id = 68 * [None] + allocation_network_id[0] = 1 + allocation_network_id[65] = 3 + allocation_network_id[66] = 5 + allocation_network_id[67] = 7 + + lines = node.geometry_from_connectivity(from_id.tolist(), to_id.tolist()) + edge = ribasim.Edge( + df=gpd.GeoDataFrame( + data={ + "from_node_id": from_id, + "to_node_id": to_id, + "edge_type": edge_type, + "allocation_network_id": allocation_network_id, + }, + geometry=lines, + crs="EPSG:28992", + ) + ) + + # Setup the basins: + profile = pd.DataFrame( + data={ + "node_id": [ + 2, + 2, + 4, + 4, + 6, + 6, + 8, + 8, + 10, + 10, + 12, + 12, + 16, + 16, + 18, + 18, + 25, + 25, + 28, + 28, + 31, + 31, + 35, + 35, + 40, + 40, + 42, + 42, + 44, + 44, + 48, + 48, + 50, + 50, + 54, + 54, + ], + "area": [ + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 100000.0, + 100000.0, + 100000.0, + 100000.0, + 100000.0, + 100000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + 1000.0, + ], + "level": [ + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + 0.0, + 1.0, + ], + } + ) + + static = pd.DataFrame( + data={ + "node_id": [ + 2, + 4, + 6, + 8, + 10, + 12, + 16, + 18, + 25, + 28, + 31, + 35, + 40, + 44, + 42, + 48, + 50, + 54, + ], + "drainage": 0.0, + "potential_evaporation": 0.0, + "infiltration": 0.0, + "precipitation": 0.0, + "urban_runoff": 0.0, + } + ) + + state = pd.DataFrame( + data={ + "node_id": [ + 2, + 4, + 6, + 8, + 10, + 12, + 16, + 18, + 25, + 28, + 31, + 35, + 40, + 42, + 44, + 48, + 50, + 54, + ], + "level": [ + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 10.0, + 10.0, + 10.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + 1.0, + ], + } + ) + + basin = ribasim.Basin( + profile=profile, + static=static, + state=state, + ) + + # Setup the discrete control: + condition = pd.DataFrame( + data={ + "node_id": [33], + "listen_feature_id": [25], + "variable": ["level"], + "greater_than": [0.003], + } + ) + + logic = pd.DataFrame( + data={ + "node_id": [33, 33], + "truth_state": ["F", "T"], + "control_state": ["A", "B"], + } + ) + + discrete_control = ribasim.DiscreteControl(condition=condition, logic=logic) + + # Setup flow boundary + flow_boundary = ribasim.FlowBoundary( + static=pd.DataFrame(data={"node_id": [1], "flow_rate": [1.0]}) + ) + + # Setup fractional flow + fractional_flow = ribasim.FractionalFlow( + static=pd.DataFrame( + data={ + "node_id": [27, 30, 27, 30], + "fraction": [0.25, 0.75, 0.75, 0.25], + "control_state": ["A", "A", "B", "B"], + } + ) + ) + + # Setup linear resistance + linear_resistance = ribasim.LinearResistance( + static=pd.DataFrame(data={"node_id": [3, 5, 7, 9], "resistance": 0.001}) + ) + + # Setup outlet + outlet = ribasim.Outlet( + static=pd.DataFrame( + data={ + "node_id": [13, 17, 23, 36, 41, 43, 55], + "flow_rate": [3.0, 3.0, 3.0, 0.003, 0.003, 0.003, 0.003], + "max_flow_rate": 3.0, + } + ) + ) + + # Setup pump + pump = ribasim.Pump( + static=pd.DataFrame( + data={ + "node_id": [15, 39, 49, 11, 24, 38], + "flow_rate": [4.0e00, 4.0e-03, 4.0e-03, 1.0e-03, 1.0e-03, 1.0e-03], + "max_flow_rate": [4.0, 0.004, 0.004, 1.0, 1.0, 1.0], + } + ) + ) + + # Setup tabulated rating curve + rating_curve = ribasim.TabulatedRatingCurve( + static=pd.DataFrame( + data={ + "node_id": [26, 26, 46, 46, 47, 47, 52, 52], + "level": [0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0], + "flow_rate": [ + 0.0e00, + 1.0e-04, + 0.0e00, + 2.0e00, + 0.0e00, + 2.0e00, + 0.0e00, + 2.0e00, + ], + } + ) + ) + + # Setup terminal node + terminal = ribasim.Terminal(static=pd.DataFrame(data={"node_id": [14, 19, 37, 56]})) + + # Setup the user + user = ribasim.User( + static=pd.DataFrame( + data={ + "node_id": [20, 21, 22, 29, 34, 45, 51, 53, 57], + "demand": [ + 4.0e00, + 5.0e00, + 3.0e00, + 1.0e-03, + 1.0e-03, + 1.0e-03, + 1.0e-03, + 1.0e-03, + 1.0e-03, + ], + "return_factor": 0.9, + "min_level": 0.9, + "priority": [2, 1, 2, 1, 2, 1, 3, 3, 2], + } + ), + time=pd.DataFrame( + data={ + "node_id": [32, 32], + "time": ["2020-01-01 00:00:00", "2021-01-01 00:00:00"], + "demand": [0.001, 0.002], + "return_factor": 0.9, + "min_level": 0.9, + "priority": 1, + } + ), + ) + + # Setup allocation: + allocation = ribasim.Allocation(use_allocation=True, timestep=86400) + + model = ribasim.Model( + network=ribasim.Network(node=node, edge=edge), + basin=basin, + discrete_control=discrete_control, + flow_boundary=flow_boundary, + fractional_flow=fractional_flow, + linear_resistance=linear_resistance, + outlet=outlet, + pump=pump, + terminal=terminal, + user=user, + tabulated_rating_curve=rating_curve, + allocation=allocation, + starttime="2020-01-01 00:00:00", + endtime="2021-01-01 00:00:00", + ) + + return model From 21d54c05d69bf4987ebce1d055c6bff822c8c105 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Fri, 2 Feb 2024 13:55:30 +0100 Subject: [PATCH 4/4] Remove confusing sentence for PidControl (#1009) People are reading this as if the outneighbor of a PidControl can only be a Pump, but that is not the case. Since it's quite a detailed and perhaps obvious point, this removes the extra line. --- docs/core/usage.qmd | 2 -- 1 file changed, 2 deletions(-) diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index 6c7e3a529..92a9c64fa 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -655,8 +655,6 @@ control_state | String The PidControl node controls the level in a basin by continuously controlling the flow rate of a connected pump or outlet. See also [PID controller](https://en.wikipedia.org/wiki/PID_controller). When A PidControl node is made inactive, the node under its control retains the last flow rate value, and the error integral is reset to 0. -The pump must be an outneighbor of the controlled basin, the outlet must be an inneighbor. Thus for the same coefficient values in `proportional`, `integral` and `derivative`, the pump and outlet have the oppositve effect on the controlled basin. - In the future controlling the flow on a particular edge could be supported. column | type | unit | restriction