diff --git a/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_Ribasim_BuildLibribasimLinux.xml b/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_Ribasim_BuildLibribasimLinux.xml index edba28fe3..d95d199aa 100644 --- a/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_Ribasim_BuildLibribasimLinux.xml +++ b/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_Ribasim_BuildLibribasimLinux.xml @@ -1,5 +1,5 @@ - + Build libribasim - Linux diff --git a/Manifest.toml b/Manifest.toml index a85784c93..55616f291 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.0-rc2" manifest_format = "2.0" -project_hash = "ffe3d5f3606c1f85a4bca07ee2df86a28560d8a2" +project_hash = "e92257faf6a6264eb6078cd603308efd2e42b380" [[deps.ADTypes]] git-tree-sha1 = "332e5d7baeff8497b923b730b994fa480601efc7" @@ -1050,9 +1050,9 @@ version = "6.59.3" [[deps.PackageCompiler]] deps = ["Artifacts", "Glob", "LazyArtifacts", "Libdl", "Pkg", "Printf", "RelocatableFolders", "TOML", "UUIDs", "p7zip_jll"] -git-tree-sha1 = "f9392ab72832f4315220a853747ff3dba758c9d1" +git-tree-sha1 = "8b880733c61c4a99a069302390de2d057ee166ce" uuid = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" -version = "2.1.15" +version = "2.1.16" [[deps.PackageExtensionCompat]] git-tree-sha1 = "fb28e33b8a95c4cee25ce296c817d89cc2e53518" diff --git a/core/src/allocation.jl b/core/src/allocation.jl index 72551dc82..ac6a0745f 100644 --- a/core/src/allocation.jl +++ b/core/src/allocation.jl @@ -3,17 +3,17 @@ 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) = p + (; graph, basin, fractional_flow) = p node_ids = graph[].node_ids[allocation_network_id] used_nodes = Set{NodeID}() - for node_id in node_ids + 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] + if node_type in [:user, :basin, :terminal] push!(used_nodes, node_id) - elseif count(x -> true, inoutflow_ids(graph, node_id)) > 2 - # use count since the length of the iterator is unknown + elseif has_fractional_flow_outneighbors push!(used_nodes, node_id) end end @@ -378,7 +378,7 @@ function add_constraints_capacity!( F = problem[:F] edge_ids = graph[].edge_ids[allocation_network_id] edge_ids_finite_capacity = Tuple{NodeID, NodeID}[] - for (i, edge) in enumerate(edge_ids) + for edge in edge_ids if !isinf(capacity[edge...]) push!(edge_ids_finite_capacity, edge) end @@ -500,13 +500,12 @@ function add_constraints_user_returnflow!( node_id for node_id in node_ids if graph[node_id].type == :user && !isempty(outflow_ids_allocation(graph, node_id)) ] - problem[:return_flow] = JuMP.@constraint( problem, [node_id_user = node_ids_user_with_returnflow], - F[Int(node_id_user), Int(only(outflow_ids_allocation(graph, node_id_user)))] <= - user.return_factor[findsorted(user.node_id, node_id)] * - F[Int(only(inflow_ids_allocation(graph, node_id_user))), Int(node_iduser)], + F[(node_id_user, only(outflow_ids_allocation(graph, node_id_user)))] <= + user.return_factor[findsorted(user.node_id, node_id_user)] * + F[(only(inflow_ids_allocation(graph, node_id_user)), node_id_user)], base_name = "return_flow", ) return nothing @@ -577,6 +576,57 @@ function add_constraints_absolute_value!( return nothing end +""" +Add the fractional flow constraints to the allocation problem. +The constraint indices are allocation edges over a fractional flow node. + +Constraint: +flow after fractional_flow node <= fraction * inflow +""" +function add_constraints_fractional_flow!( + problem::JuMP.Model, + p::Parameters, + allocation_network_id::Int, +)::Nothing + (; graph, fractional_flow) = p + F = problem[:F] + node_ids = graph[].node_ids[allocation_network_id] + + edges_to_fractional_flow = Tuple{NodeID, NodeID}[] + fractions = Dict{Tuple{NodeID, NodeID}, Float64}() + inflows = Dict{NodeID, JuMP.AffExpr}() + for node_id in node_ids + for outflow_id_ in outflow_ids(graph, node_id) + if graph[outflow_id_].type == :fractional_flow + # The fractional flow nodes themselves are not represented in + # the allocation graph + dst_id = outflow_id(graph, outflow_id_) + # For now only consider fractional flow nodes which end in a basin + if haskey(graph, node_id, dst_id) && graph[dst_id].type == :basin + edge = (node_id, dst_id) + push!(edges_to_fractional_flow, edge) + node_idx = findsorted(fractional_flow.node_id, outflow_id_) + fractions[edge] = fractional_flow.fraction[node_idx] + inflows[node_id] = sum([ + F[(inflow_id_, node_id)] for + inflow_id_ in inflow_ids(graph, node_id) + ]) + end + end + end + end + + if !isempty(edges_to_fractional_flow) + problem[:fractional_flow] = JuMP.@constraint( + problem, + [edge = edges_to_fractional_flow], + F[edge] <= fractions[edge] * inflows[edge[1]], + base_name = "fractional_flow" + ) + end + return nothing +end + """ Construct the allocation problem for the current subnetwork as a JuMP.jl model. """ @@ -600,7 +650,7 @@ function allocation_problem( add_constraints_flow_conservation!(problem, p, allocation_network_id) add_constraints_user_returnflow!(problem, p, allocation_network_id) add_constraints_absolute_value!(problem, p, allocation_network_id, config) - # TODO: The fractional flow constraints + add_constraints_fractional_flow!(problem, p, allocation_network_id) return problem end diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 5da20caa1..d5d5d9ba3 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -32,10 +32,6 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model try parameters = Parameters(db, config) - if !valid_n_neighbors(parameters) - error("Invalid number of connections for certain node types.") - end - if !valid_discrete_control(parameters, config) error("Invalid discrete control state definition(s).") end @@ -485,6 +481,51 @@ function discrete_control_affect!( return control_state_change end +function get_allocation_model( + 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 + end + return nothing +end + +""" +Update the fractional flow fractions in an allocation problem. +""" +function set_fractional_flow_in_allocation!( + p::Parameters, + node_id::NodeID, + fraction::Number, +)::Nothing + (; graph) = p + + allocation_network_id = graph[node_id].allocation_network_id + # Get the allocation model this fractional flow node is in + allocation_model = get_allocation_model(p, allocation_network_id) + if !isnothing(allocation_model) + problem = allocation_model.problem + # The allocation edge which jumps over the fractional flow node + edge = (inflow_id(graph, node_id), outflow_id(graph, node_id)) + if haskey(graph, edge...) + # The constraint for this fractional flow node + if edge in keys(problem[:fractional_flow]) + constraint = problem[:fractional_flow][edge] + + # Set the new fraction on all inflow terms in the constraint + for inflow_id in inflow_ids_allocation(graph, edge[1]) + flow = problem[:F][(inflow_id, edge[1])] + JuMP.set_normalized_coefficient(constraint, flow, -fraction) + end + end + end + end + return nothing +end + function set_control_params!(p::Parameters, node_id::NodeID, control_state::String) node = getfield(p, p.graph[node_id].type) idx = searchsortedfirst(node.node_id, node_id) @@ -495,6 +536,11 @@ function set_control_params!(p::Parameters, node_id::NodeID, control_state::Stri vec = get_tmp(getfield(node, field), 0) vec[idx] = value end + + # Set new fractional flow fractions in allocation problem + if node isa FractionalFlow && field == :fraction + set_fractional_flow_in_allocation!(p, node_id, value) + end end end diff --git a/core/src/create.jl b/core/src/create.jl index aef043663..46a5f9ada 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -846,6 +846,11 @@ function Parameters(db::DB, config::Config)::Parameters Dict{Int, Symbol}(), subgrid_level, ) + + if !valid_n_neighbors(p) + error("Invalid number of connections for certain node types.") + end + # Allocation data structures if config.allocation.use_allocation generate_allocation_models!(p, config) diff --git a/core/src/utils.jl b/core/src/utils.jl index cfa1d254d..5040f786c 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -1206,7 +1206,7 @@ is_flow_constraining(node::AbstractParameterNode) = hasfield(typeof(node), :max_ """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]) + (nameof(typeof(node)) ∈ [:Pump, :Outlet, :TabulatedRatingCurve, :FractionalFlow]) """Find out whether a path exists between a start node and end node in the given allocation graph.""" function allocation_path_exists_in_graph( diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 12b756464..ff48730e6 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -109,3 +109,59 @@ end @test objective.terms[F[(NodeID(4), NodeID(5))]] ≈ 62.585499316005475 @test objective.terms[F[(NodeID(2), NodeID(4))]] ≈ 62.585499316005475 end + +@testitem "Allocation with controlled fractional flow" begin + using DataFrames + using Ribasim: NodeID + using OrdinaryDiffEq: solve! + using JuMP + + toml_path = normpath( + @__DIR__, + "../../generated_testmodels/fractional_flow_subnetwork/ribasim.toml", + ) + model = Ribasim.BMI.initialize(Ribasim.Model, toml_path) + problem = model.integrator.p.allocation_models[1].problem + F = problem[:F] + @test JuMP.normalized_coefficient( + problem[:fractional_flow][(NodeID(3), NodeID(5))], + F[(NodeID(2), NodeID(3))], + ) ≈ -0.75 + @test JuMP.normalized_coefficient( + problem[:fractional_flow][(NodeID(3), NodeID(8))], + F[(NodeID(2), NodeID(3))], + ) ≈ -0.25 + + solve!(model) + record_allocation = DataFrame(model.integrator.p.user.record) + record_control = model.integrator.p.discrete_control.record + groups = groupby(record_allocation, [:user_node_id, :priority]) + fractional_flow = model.integrator.p.fractional_flow + (; control_mapping) = fractional_flow + t_control = record_control.time[2] + + allocated_6_before = groups[(6, 1)][groups[(6, 1)].time .< t_control, :].allocated + allocated_9_before = groups[(9, 1)][groups[(9, 1)].time .< t_control, :].allocated + allocated_6_after = groups[(6, 1)][groups[(6, 1)].time .> t_control, :].allocated + allocated_9_after = groups[(9, 1)][groups[(9, 1)].time .> t_control, :].allocated + @test all( + allocated_9_before ./ allocated_6_before .<= + control_mapping[(NodeID(7), "A")].fraction / + control_mapping[(NodeID(4), "A")].fraction, + ) + @test all(allocated_9_after ./ allocated_6_after .<= 1.0) + + @test record_control.truth_state == ["F", "T"] + @test record_control.control_state == ["A", "B"] + + fractional_flow_constraints = + model.integrator.p.allocation_models[1].problem[:fractional_flow] + @test JuMP.normalized_coefficient( + problem[:fractional_flow][(NodeID(3), NodeID(5))], + F[(NodeID(2), NodeID(3))], + ) ≈ -0.75 + @test JuMP.normalized_coefficient( + problem[:fractional_flow][(NodeID(3), NodeID(8))], + F[(NodeID(2), NodeID(3))], + ) ≈ -0.25 +end diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index 0a03bb112..15de4dcca 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -169,6 +169,17 @@ end Sys.isapple() end +@testitem "allocation example model" begin + using SciMLBase: successful_retcode + + toml_path = + normpath(@__DIR__, "../../generated_testmodels/allocation_example/ribasim.toml") + @test ispath(toml_path) + model = Ribasim.run(Ribasim.Config(toml_path)) + @test model isa Ribasim.Model + @test successful_retcode(model) +end + @testitem "sparse and AD/FDM jac solver options" begin using SciMLBase: successful_retcode diff --git a/docs/Project.toml b/docs/Project.toml index 1af8db056..a15dc604d 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -22,12 +22,15 @@ Ribasim = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635" [compat] Configurations = "0.17" +DataFrames = "1" Dates = "<0.0.1,1" Documenter = "0.27,1" DocumenterMarkdown = "0.2" +IJulia = "1" InteractiveUtils = "<0.0.1,1" JSON3 = "1.12" Legolas = "0.5" Logging = "<0.0.1,1" +MarkdownTables = "1" OrderedCollections = "1.6" julia = "1.10" diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index 744589add..f80c79a9c 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -166,10 +166,6 @@ $$ {#eq-fractionalflowconstraint} - Flow sign: Furthermore there are the non-negativity constraints for the flows and allocations, see [The optimization variables](allocation.qmd#the-optimization-variables). -:::{.callout-note} -Currently the fractional flow constraints are not taken into account in the implementation. -::: - ## Final notes on the allocation problem ### Users using their own return flow diff --git a/docs/core/equations_files/validation_tables/connectivity.jl b/docs/core/equations_files/validation_tables/connectivity.jl deleted file mode 100644 index e69de29bb..000000000 diff --git a/docs/python/examples.ipynb b/docs/python/examples.ipynb index 40f863117..4e00ab232 100644 --- a/docs/python/examples.ipynb +++ b/docs/python/examples.ipynb @@ -1425,6 +1425,504 @@ "ax.plot(times, target_levels, color=\"k\", ls=\":\", label=\"target level\")\n", "pass" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Model with allocation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the nodes:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "xy = np.array(\n", + " [\n", + " (0.0, 0.0), # 1: FlowBoundary\n", + " (1.0, 0.0), # 2: Basin\n", + " (1.0, 1.0), # 3: User\n", + " (2.0, 0.0), # 4: LinearResistance\n", + " (3.0, 0.0), # 5: Basin\n", + " (3.0, 1.0), # 6: User\n", + " (4.0, 0.0), # 7: TabulatedRatingCurve\n", + " (4.5, 0.0), # 8: FractionalFlow\n", + " (4.5, 0.5), # 9: FractionalFlow\n", + " (5.0, 0.0), # 10: Terminal\n", + " (4.5, 0.25), # 11: DiscreteControl\n", + " (4.5, 1.0), # 12: Basin\n", + " (5.0, 1.0), # 13: User\n", + " ]\n", + ")\n", + "node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n", + "\n", + "node_type = [\n", + " \"FlowBoundary\",\n", + " \"Basin\",\n", + " \"User\",\n", + " \"LinearResistance\",\n", + " \"Basin\",\n", + " \"User\",\n", + " \"TabulatedRatingCurve\",\n", + " \"FractionalFlow\",\n", + " \"FractionalFlow\",\n", + " \"Terminal\",\n", + " \"DiscreteControl\",\n", + " \"Basin\",\n", + " \"User\",\n", + "]\n", + "\n", + "# All nodes belong to allocation network id 1\n", + "node = ribasim.Node(\n", + " df=gpd.GeoDataFrame(\n", + " data={\"type\": node_type, \"allocation_network_id\": 1},\n", + " index=pd.Index(np.arange(len(xy)) + 1, name=\"fid\"),\n", + " geometry=node_xy,\n", + " crs=\"EPSG:28992\",\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the edges:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from_id = np.array(\n", + " [1, 2, 2, 4, 5, 5, 7, 3, 6, 7, 8, 9, 12, 13, 11, 11],\n", + " dtype=np.int64,\n", + ")\n", + "to_id = np.array(\n", + " [2, 3, 4, 5, 6, 7, 8, 2, 5, 9, 10, 12, 13, 10, 8, 9],\n", + " dtype=np.int64,\n", + ")\n", + "# Denote the first edge, 1 => 2, as a source edge for\n", + "# allocation network 1\n", + "allocation_network_id = len(from_id) * [None]\n", + "allocation_network_id[0] = 1\n", + "lines = node.geometry_from_connectivity(from_id, to_id)\n", + "edge = ribasim.Edge(\n", + " df=gpd.GeoDataFrame(\n", + " data={\n", + " \"from_node_id\": from_id,\n", + " \"to_node_id\": to_id,\n", + " \"edge_type\": (len(from_id) - 2) * [\"flow\"] + 2 * [\"control\"],\n", + " \"allocation_network_id\": allocation_network_id,\n", + " },\n", + " geometry=lines,\n", + " crs=\"EPSG:28992\",\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the basins:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "profile = pd.DataFrame(\n", + " data={\n", + " \"node_id\": [2, 2, 5, 5, 12, 12],\n", + " \"area\": 300_000.0,\n", + " \"level\": 3 * [0.0, 1.0],\n", + " }\n", + ")\n", + "\n", + "static = pd.DataFrame(\n", + " data={\n", + " \"node_id\": [2, 5, 12],\n", + " \"drainage\": 0.0,\n", + " \"potential_evaporation\": 0.0,\n", + " \"infiltration\": 0.0,\n", + " \"precipitation\": 0.0,\n", + " \"urban_runoff\": 0.0,\n", + " }\n", + ")\n", + "\n", + "state = pd.DataFrame(data={\"node_id\": [2, 5, 12], \"level\": 1.0})\n", + "\n", + "basin = ribasim.Basin(profile=profile, static=static, state=state)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the flow boundary:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "flow_boundary = ribasim.FlowBoundary(\n", + " static=pd.DataFrame(\n", + " data={\n", + " \"node_id\": [1],\n", + " \"flow_rate\": 2.0,\n", + " }\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the linear resistance:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "linear_resistance = ribasim.LinearResistance(\n", + " static=pd.DataFrame(\n", + " data={\n", + " \"node_id\": [4],\n", + " \"resistance\": 0.06,\n", + " }\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the tabulated rating curve:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tabulated_rating_curve = ribasim.TabulatedRatingCurve(\n", + " static=pd.DataFrame(\n", + " data={\n", + " \"node_id\": 7,\n", + " \"level\": [0.0, 0.5, 1.0],\n", + " \"discharge\": [0.0, 0.0, 2.0],\n", + " }\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the fractional flow:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fractional_flow = ribasim.FractionalFlow(\n", + " static=pd.DataFrame(\n", + " data={\n", + " \"node_id\": [8, 8, 9, 9],\n", + " \"fraction\": [0.6, 0.9, 0.4, 0.1],\n", + " \"control_state\": [\"divert\", \"close\", \"divert\", \"close\"],\n", + " }\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the terminal:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terminal = ribasim.Terminal(\n", + " static=pd.DataFrame(\n", + " data={\n", + " \"node_id\": [10],\n", + " }\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the discrete control:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "condition = pd.DataFrame(\n", + " data={\n", + " \"node_id\": [11],\n", + " \"listen_feature_id\": 5,\n", + " \"variable\": \"level\",\n", + " \"greater_than\": 0.52,\n", + " }\n", + ")\n", + "\n", + "logic = pd.DataFrame(\n", + " data={\n", + " \"node_id\": 11,\n", + " \"truth_state\": [\"T\", \"F\"],\n", + " \"control_state\": [\"divert\", \"close\"],\n", + " }\n", + ")\n", + "\n", + "discrete_control = ribasim.DiscreteControl(condition=condition, logic=logic)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the users:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "user = ribasim.User(\n", + " static=pd.DataFrame(\n", + " data={\n", + " \"node_id\": [6, 13],\n", + " \"demand\": [1.5, 1.0],\n", + " \"return_factor\": 0.0,\n", + " \"min_level\": -1.0,\n", + " \"priority\": [1, 3],\n", + " }\n", + " ),\n", + " time=pd.DataFrame(\n", + " data={\n", + " \"node_id\": [3, 3, 3, 3],\n", + " \"demand\": [0.0, 1.0, 1.2, 1.2],\n", + " \"priority\": [1, 1, 2, 2],\n", + " \"return_factor\": 0.0,\n", + " \"min_level\": -1.0,\n", + " \"time\": 2 * [\"2020-01-01 00:00:00\", \"2020-01-20 00:00:00\"],\n", + " }\n", + " ),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the allocation:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "allocation = ribasim.Allocation(use_allocation=True, timestep=86400)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup a model:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model = ribasim.Model(\n", + " network=ribasim.Network(\n", + " node=node,\n", + " edge=edge,\n", + " ),\n", + " basin=basin,\n", + " flow_boundary=flow_boundary,\n", + " linear_resistance=linear_resistance,\n", + " tabulated_rating_curve=tabulated_rating_curve,\n", + " terminal=terminal,\n", + " user=user,\n", + " discrete_control=discrete_control,\n", + " fractional_flow=fractional_flow,\n", + " allocation=allocation,\n", + " starttime=\"2020-01-01 00:00:00\",\n", + " endtime=\"2020-01-20 00:00:00\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's take a look at the model:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Write the model to a TOML and GeoPackage:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "datadir = Path(\"data\")\n", + "model.write(datadir / \"allocation_example/ribasim.toml\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# | include: false\n", + "from subprocess import run\n", + "\n", + "run(\n", + " [\n", + " \"julia\",\n", + " \"--project=../../core\",\n", + " \"--eval\",\n", + " f'using Ribasim; Ribasim.run(\"{datadir.as_posix()}/allocation_example/ribasim.toml\")',\n", + " ],\n", + " check=True,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now run the model with `ribasim allocation_example/ribasim.toml`.\n", + "After running the model, read back the results:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.ticker as plticker\n", + "\n", + "df_allocation = pd.read_feather(datadir / \"allocation_example/results/allocation.arrow\")\n", + "df_allocation_wide = df_allocation.pivot_table(\n", + " index=\"time\",\n", + " columns=[\"user_node_id\", \"priority\"],\n", + " values=[\"demand\", \"allocated\", \"abstracted\"],\n", + ")\n", + "df_allocation_wide = df_allocation_wide.loc[:, (df_allocation_wide != 0).any(axis=0)]\n", + "\n", + "fig, axs = plt.subplots(1, 3, figsize=(8, 5))\n", + "fig.suptitle(f\"Objective type: {model.allocation.objective_type}\", fontweight=\"bold\")\n", + "\n", + "df_allocation_wide[\"demand\"].plot(ax=axs[0], ls=\":\")\n", + "df_allocation_wide[\"allocated\"].plot(ax=axs[1], ls=\"--\")\n", + "df_allocation_wide[\"abstracted\"].plot(ax=axs[2])\n", + "\n", + "fig.tight_layout()\n", + "loc = plticker.MultipleLocator(2)\n", + "\n", + "axs[0].set_ylabel(\"level [m]\")\n", + "\n", + "for ax, title in zip(axs, [\"Demand\", \"Allocated\", \"Abstracted\"]):\n", + " ax.set_title(title)\n", + " ax.set_ylim(0.0, 1.6)\n", + " ax.xaxis.set_major_locator(loc)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Some things to note about this plot:\n", + "\n", + "- Abstraction behaves somewhat erratically at the start of the simulation. This is because allocation is based on flows computed in the physical layer, and at the start of the simulation these are not known yet.\n", + "\n", + "- Although there is a plotted line for abstraction per priority, abstraction is actually accumulated over all priorities per user." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df_basin = pd.read_feather(datadir / \"allocation_example/results/basin.arrow\")\n", + "df_basin_wide = df_basin.pivot_table(\n", + " index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n", + ")\n", + "\n", + "ax = df_basin_wide[\"level\"].plot()\n", + "ax.set_title(\"Basin levels\")\n", + "ax.set_ylabel(\"level [m]\")" + ] } ], "metadata": { diff --git a/pixi.toml b/pixi.toml index f6118107e..e5e01c9e1 100644 --- a/pixi.toml +++ b/pixi.toml @@ -37,7 +37,7 @@ quarto-preview = { cmd = "export QUARTO_PYTHON=python && quarto preview docs", d "quartodoc-build", ] } quarto-check = { cmd = "quarto check all", depends_on = ["quartodoc-build"] } -quarto-render = { cmd = "export QUARTO_PYTHON=python && quarto render docs --to html --execute", depends_on = [ +quarto-render = { cmd = "julia --project=docs --eval 'using Pkg; Pkg.build(\"IJulia\")' && export QUARTO_PYTHON=python && quarto render docs --to html --execute", depends_on = [ "quartodoc-build", ] } docs = { depends_on = ["build-julia-docs", "quarto-preview"] } diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index cbc8b295f..7de7eb347 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -6,8 +6,10 @@ import ribasim_testmodels from ribasim_testmodels.allocation import ( - minimal_subnetwork_model, + allocation_example_model, # looped_subnetwork_model, + fractional_flow_subnetwork_model, + minimal_subnetwork_model, subnetwork_model, user_model, ) @@ -49,6 +51,7 @@ from ribasim_testmodels.trivial import trivial_model __all__ = [ + "allocation_example_model", "backwater_model", "basic_model", "basic_arrow_model", @@ -78,6 +81,7 @@ "user_model", "subnetwork_model", "minimal_subnetwork_model", + "fractional_flow_subnetwork_model", # Disable until this issue is resolved: # https://github.com/Deltares/Ribasim/issues/692 # "looped_subnetwork_model", diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 33606be89..c3416c4d4 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -667,3 +667,401 @@ def minimal_subnetwork_model(): ) return model + + +def fractional_flow_subnetwork_model(): + """Create a small subnetwork that contains fractional flow nodes.""" + + xy = np.array( + [ + (0.0, 0.0), # 1: FlowBoundary + (0.0, 1.0), # 2: Basin + (0.0, 2.0), # 3: TabulatedRatingCurve + (-1.0, 3.0), # 4: FractionalFlow + (-2.0, 4.0), # 5: Basin + (-3.0, 5.0), # 6: User + (1.0, 3.0), # 7: FractionalFlow + (2.0, 4.0), # 8: Basin + (3.0, 5.0), # 9: User + (-1.0, 2.0), # 10: DiscreteControl + ] + ) + node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1]) + + node_type = [ + "FlowBoundary", + "Basin", + "TabulatedRatingCurve", + "FractionalFlow", + "Basin", + "User", + "FractionalFlow", + "Basin", + "User", + "DiscreteControl", + ] + + # 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}, + 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, 3, 7, 8, 9, 10, 10], + dtype=np.int64, + ) + to_id = np.array( + [2, 3, 4, 5, 6, 5, 7, 8, 9, 8, 4, 7], + dtype=np.int64, + ) + allocation_network_id = len(from_id) * [None] + allocation_network_id[0] = 1 + lines = node.geometry_from_connectivity(from_id, to_id) + edge = ribasim.Edge( + df=gpd.GeoDataFrame( + data={ + "from_node_id": from_id, + "to_node_id": to_id, + "edge_type": (len(from_id) - 2) * ["flow"] + 2 * ["control"], + "allocation_network_id": allocation_network_id, + }, + geometry=lines, + crs="EPSG:28992", + ) + ) + + # Setup the basins: + profile = pd.DataFrame( + data={ + "node_id": [2, 2, 5, 5, 8, 8], + "area": 1000.0, + "level": 3 * [0.0, 1.0], + } + ) + + static = pd.DataFrame( + data={ + "node_id": [2, 5, 8], + "drainage": 0.0, + "potential_evaporation": 0.0, + "infiltration": 0.0, + "precipitation": 0.0, + "urban_runoff": 0.0, + } + ) + + state = pd.DataFrame(data={"node_id": [2, 5, 8], "level": 1.0}) + + basin = ribasim.Basin(profile=profile, static=static, state=state) + + # Setup the flow boundary: + flow_boundary = ribasim.FlowBoundary( + time=pd.DataFrame( + data={ + "node_id": [1, 1], + "flow_rate": [2.0e-3, 4.0e-3], + "time": ["2020-01-01 00:00:00", "2021-01-01 00:00:00"], + } + ) + ) + + # Setup the tabulated rating curve: + rating_curve = ribasim.TabulatedRatingCurve( + static=pd.DataFrame( + data={ + "node_id": [3, 3], + "level": [0.0, 1.0], + "discharge": [0.0, 1e-4], + } + ) + ) + + # Setup the users: + user = ribasim.User( + static=pd.DataFrame( + data={ + "node_id": [6], + "demand": 1.0e-3, + "return_factor": 0.9, + "min_level": 0.9, + "priority": 1, + } + ), + time=pd.DataFrame( + data={ + "time": ["2020-01-01 00:00:00", "2021-01-01 00:00:00"], + "node_id": 9, + "demand": [1.0e-3, 2.0e-3], + "return_factor": 0.9, + "min_level": 0.9, + "priority": 1, + } + ), + ) + + # Setup allocation: + allocation = ribasim.Allocation(use_allocation=True, timestep=86400) + + # Setup fractional flows: + fractional_flow = ribasim.FractionalFlow( + static=pd.DataFrame( + data={ + "node_id": [4, 7, 4, 7], + "fraction": [0.25, 0.75, 0.75, 0.25], + "control_state": ["A", "A", "B", "B"], + } + ) + ) + + # Setup discrete control: + condition = pd.DataFrame( + data={ + "node_id": [10], + "listen_feature_id": [1], + "variable": "flow_rate", + "greater_than": [3.0e-3], + } + ) + + logic = pd.DataFrame( + data={ + "node_id": [10, 10], + "truth_state": ["F", "T"], + "control_state": ["A", "B"], + } + ) + + discrete_control = ribasim.DiscreteControl(condition=condition, logic=logic) + + model = ribasim.Model( + network=ribasim.Network( + node=node, + edge=edge, + ), + basin=basin, + flow_boundary=flow_boundary, + tabulated_rating_curve=rating_curve, + user=user, + allocation=allocation, + fractional_flow=fractional_flow, + discrete_control=discrete_control, + starttime="2020-01-01 00:00:00", + endtime="2021-01-01 00:00:00", + ) + + return model + + +def allocation_example_model(): + """Generate a model that is used as an example of allocation in the docs.""" + + xy = np.array( + [ + (0.0, 0.0), # 1: FlowBoundary + (1.0, 0.0), # 2: Basin + (1.0, 1.0), # 3: User + (2.0, 0.0), # 4: LinearResistance + (3.0, 0.0), # 5: Basin + (3.0, 1.0), # 6: User + (4.0, 0.0), # 7: TabulatedRatingCurve + (4.5, 0.0), # 8: FractionalFlow + (4.5, 0.5), # 9: FractionalFlow + (5.0, 0.0), # 10: Terminal + (4.5, 0.25), # 11: DiscreteControl + (4.5, 1.0), # 12: Basin + (5.0, 1.0), # 13: User + ] + ) + node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1]) + + node_type = [ + "FlowBoundary", + "Basin", + "User", + "LinearResistance", + "Basin", + "User", + "TabulatedRatingCurve", + "FractionalFlow", + "FractionalFlow", + "Terminal", + "DiscreteControl", + "Basin", + "User", + ] + + # All nodes belong to allocation network id 1 + node = ribasim.Node( + df=gpd.GeoDataFrame( + data={"type": node_type, "allocation_network_id": 1}, + index=pd.Index(np.arange(len(xy)) + 1, name="fid"), + geometry=node_xy, + crs="EPSG:28992", + ) + ) + + from_id = np.array( + [1, 2, 2, 4, 5, 5, 7, 3, 6, 7, 8, 9, 12, 13, 11, 11], + dtype=np.int64, + ) + to_id = np.array( + [2, 3, 4, 5, 6, 7, 8, 2, 5, 9, 10, 12, 13, 10, 8, 9], + dtype=np.int64, + ) + # 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 + lines = node.geometry_from_connectivity(from_id, to_id) + edge = ribasim.Edge( + df=gpd.GeoDataFrame( + data={ + "from_node_id": from_id, + "to_node_id": to_id, + "edge_type": (len(from_id) - 2) * ["flow"] + 2 * ["control"], + "allocation_network_id": allocation_network_id, + }, + geometry=lines, + crs="EPSG:28992", + ) + ) + + # Setup the basins: + profile = pd.DataFrame( + data={ + "node_id": [2, 2, 5, 5, 12, 12], + "area": 300_000.0, + "level": 3 * [0.0, 1.0], + } + ) + + static = pd.DataFrame( + data={ + "node_id": [2, 5, 12], + "drainage": 0.0, + "potential_evaporation": 0.0, + "infiltration": 0.0, + "precipitation": 0.0, + "urban_runoff": 0.0, + } + ) + + state = pd.DataFrame(data={"node_id": [2, 5, 12], "level": 1.0}) + + basin = ribasim.Basin(profile=profile, static=static, state=state) + + flow_boundary = ribasim.FlowBoundary( + static=pd.DataFrame( + data={ + "node_id": [1], + "flow_rate": 2.0, + } + ) + ) + + linear_resistance = ribasim.LinearResistance( + static=pd.DataFrame( + data={ + "node_id": [4], + "resistance": 0.06, + } + ) + ) + + tabulated_rating_curve = ribasim.TabulatedRatingCurve( + static=pd.DataFrame( + data={ + "node_id": 7, + "level": [0.0, 0.5, 1.0], + "discharge": [0.0, 0.0, 2.0], + } + ) + ) + + fractional_flow = ribasim.FractionalFlow( + static=pd.DataFrame( + data={ + "node_id": [8, 8, 9, 9], + "fraction": [0.6, 0.9, 0.4, 0.1], + "control_state": ["divert", "close", "divert", "close"], + } + ) + ) + + terminal = ribasim.Terminal( + static=pd.DataFrame( + data={ + "node_id": [10], + } + ) + ) + + condition = pd.DataFrame( + data={ + "node_id": [11], + "listen_feature_id": 5, + "variable": "level", + "greater_than": 0.52, + } + ) + + logic = pd.DataFrame( + data={ + "node_id": 11, + "truth_state": ["T", "F"], + "control_state": ["divert", "close"], + } + ) + + discrete_control = ribasim.DiscreteControl(condition=condition, logic=logic) + + user = ribasim.User( + static=pd.DataFrame( + data={ + "node_id": [6, 13], + "demand": [1.5, 1.0], + "return_factor": 0.0, + "min_level": -1.0, + "priority": [1, 3], + } + ), + time=pd.DataFrame( + data={ + "node_id": [3, 3, 3, 3], + "demand": [0.0, 1.0, 1.2, 1.2], + "priority": [1, 1, 2, 2], + "return_factor": 0.0, + "min_level": -1.0, + "time": 2 * ["2020-01-01 00:00:00", "2020-01-20 00:00:00"], + } + ), + ) + + # Setup allocation: + allocation = ribasim.Allocation(use_allocation=True, timestep=86400) + + model = ribasim.Model( + network=ribasim.Network( + node=node, + edge=edge, + ), + basin=basin, + flow_boundary=flow_boundary, + linear_resistance=linear_resistance, + tabulated_rating_curve=tabulated_rating_curve, + terminal=terminal, + user=user, + discrete_control=discrete_control, + fractional_flow=fractional_flow, + allocation=allocation, + starttime="2020-01-01 00:00:00", + endtime="2020-01-20 00:00:00", + ) + + return model