From 8888cd2ac6e636b69316dc7bc6a177536d6188d9 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Mon, 15 Apr 2024 11:28:12 +0200 Subject: [PATCH] Revert "Equal ratio allocation objective (#1366)" This reverts commit 19139f1c8cc7c539e7f18ea329b14f8226ceb79d. --- core/src/allocation_init.jl | 158 +++++++++++++++++- core/src/allocation_optim.jl | 105 ++++++++---- core/test/allocation_test.jl | 115 +++++-------- docs/core/allocation.qmd | 12 +- .../ribasim_testmodels/__init__.py | 2 - .../ribasim_testmodels/allocation.py | 97 ----------- 6 files changed, 281 insertions(+), 208 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 0f16cb2b4..d17ffc094 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -5,8 +5,7 @@ function find_subnetwork_connections!(p::Parameters)::Nothing (; 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) | - (graph[node_id, outflow_id].allocation_network_id_source == 1) + 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) @@ -406,6 +405,55 @@ function add_variables_flow_buffer!( return nothing end +""" +Certain allocation distribution types use absolute values in the objective function. +Since most optimization packages do not support the absolute value function directly, +New variables are introduced that act as the absolute value of an expression by +posing the appropriate constraints. +""" +function add_variables_absolute_value!( + problem::JuMP.Model, + p::Parameters, + allocation_network_id::Int32, +)::Nothing + (; graph, allocation) = p + (; main_network_connections) = allocation + + node_ids = graph[].node_ids[allocation_network_id] + node_ids_user_demand = NodeID[] + node_ids_level_demand = NodeID[] + node_ids_flow_demand = NodeID[] + + for node_id in node_ids + type = node_id.type + if type == NodeType.UserDemand + push!(node_ids_user_demand, node_id) + elseif type == NodeType.Basin + push!(node_ids_level_demand, node_id) + elseif has_external_demand(graph, node_id, :flow_demand)[1] + push!(node_ids_flow_demand, node_id) + end + end + + # For the main network, connections to subnetworks are treated as UserDemands + if is_main_network(allocation_network_id) + for connections_subnetwork in main_network_connections + for connection in connections_subnetwork + push!(node_ids_user_demand, connection[2]) + end + end + end + + problem[:F_abs_user_demand] = + JuMP.@variable(problem, F_abs_user_demand[node_id = node_ids_user_demand]) + problem[:F_abs_level_demand] = + JuMP.@variable(problem, F_abs_level_demand[node_id = node_ids_level_demand]) + problem[:F_abs_flow_demand] = + JuMP.@variable(problem, F_abs_flow_demand[node_id = node_ids_flow_demand]) + + return nothing +end + """ Add the flow capacity constraints to the allocation problem. Only finite capacities get a constraint. @@ -420,6 +468,7 @@ function add_constraints_capacity!( p::Parameters, allocation_network_id::Int32, )::Nothing + (; graph) = p main_network_source_edges = get_main_network_connections(p, allocation_network_id) F = problem[:F] edge_ids_finite_capacity = Tuple{NodeID, NodeID}[] @@ -643,6 +692,106 @@ function add_constraints_conservation_basin!( return nothing end +""" +Minimizing |expr| can be achieved by introducing a new variable expr_abs +and posing the following constraints: +expr_abs >= expr +expr_abs >= -expr +""" +function add_constraints_absolute_value!( + problem::JuMP.Model, + flow_per_node::Dict{NodeID, JuMP.VariableRef}, + F_abs::JuMP.Containers.DenseAxisArray, + variable_type::String, +)::Nothing + # Example demand + d = 2.0 + + node_ids = only(F_abs.axes) + + # These constraints together make sure that F_abs_* acts as the absolute + # value F_abs_* = |x| where x = F-d (here for example d = 2) + base_name = "abs_positive_$variable_type" + problem[Symbol(base_name)] = JuMP.@constraint( + problem, + [node_id = node_ids], + F_abs[node_id] >= (flow_per_node[node_id] - d), + base_name = base_name + ) + base_name = "abs_negative_$variable_type" + problem[Symbol(base_name)] = JuMP.@constraint( + problem, + [node_id = node_ids], + F_abs[node_id] >= -(flow_per_node[node_id] - d), + base_name = base_name + ) + + return nothing +end + +""" +Add constraints so that variables F_abs_user_demand act as the +absolute value of the expression comparing flow to a UserDemand to its demand. +""" +function add_constraints_absolute_value_user_demand!( + problem::JuMP.Model, + p::Parameters, +)::Nothing + (; graph) = p + + F = problem[:F] + F_abs_user_demand = problem[:F_abs_user_demand] + + flow_per_node = Dict( + node_id => F[(only(inflow_ids_allocation(graph, node_id)), node_id)] for + node_id in only(F_abs_user_demand.axes) + ) + + add_constraints_absolute_value!( + problem, + flow_per_node, + F_abs_user_demand, + "user_demand", + ) + + return nothing +end + +""" +Add constraints so that variables F_abs_level_demand act as the +absolute value of the expression comparing flow to a basin to its demand. +""" +function add_constraints_absolute_value_level_demand!(problem::JuMP.Model)::Nothing + F_basin_in = problem[:F_basin_in] + F_abs_level_demand = problem[:F_abs_level_demand] + flow_per_node = + Dict(node_id => F_basin_in[node_id] for node_id in only(F_abs_level_demand.axes)) + + add_constraints_absolute_value!(problem, flow_per_node, F_abs_level_demand, "basin") + + return nothing +end + +""" +Add constraints so that variables F_abs_flow_demand act as the +absolute value of the expression comparing flow to a flow buffer to the flow demand. +""" +function add_constraints_absolute_value_flow_demand!(problem::JuMP.Model)::Nothing + F_flow_buffer_in = problem[:F_flow_buffer_in] + F_abs_flow_demand = problem[:F_abs_flow_demand] + flow_per_node = Dict( + node_id => F_flow_buffer_in[node_id] for node_id in only(F_abs_flow_demand.axes) + ) + + add_constraints_absolute_value!( + problem, + flow_per_node, + F_abs_flow_demand, + "flow_demand", + ) + return nothing +end + """ Add the fractional flow constraints to the allocation problem. The constraint indices are allocation edges over a fractional flow node. @@ -772,6 +921,7 @@ function allocation_problem( # Add variables to problem add_variables_flow!(problem, p, allocation_network_id) add_variables_basin!(problem, p, allocation_network_id) + add_variables_absolute_value!(problem, p, allocation_network_id) add_variables_flow_buffer!(problem, p, allocation_network_id) # Add constraints to problem @@ -779,6 +929,10 @@ function allocation_problem( add_constraints_conservation_flow_demand!(problem, p, allocation_network_id) add_constraints_conservation_subnetwork!(problem, p, allocation_network_id) + add_constraints_absolute_value_user_demand!(problem, p) + add_constraints_absolute_value_flow_demand!(problem) + add_constraints_absolute_value_level_demand!(problem) + add_constraints_capacity!(problem, capacity, p, allocation_network_id) add_constraints_source!(problem, p, allocation_network_id) add_constraints_user_source!(problem, p, allocation_network_id) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 9172c13b2..b42393f75 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -1,26 +1,67 @@ @enumx OptimizationType internal_sources collect_demands allocate """ -Add an objective term `demand * (1 - flow/demand)^2`. If the absolute -value of the demand is very small, this would lead to huge coefficients, -so in that case a term of the form (flow - demand)^2 is used. +Add a term to the objective function given by the objective type, +depending in the provided flow variable and the associated demand. """ function add_objective_term!( - ex::JuMP.QuadExpr, demand::Float64, - F::JuMP.VariableRef, + constraint_abs_positive::Union{JuMP.ConstraintRef, Nothing} = nothing, + constraint_abs_negative::Union{JuMP.ConstraintRef, Nothing} = nothing, )::Nothing - if abs(demand) < 1e-5 - # Error term (F - d)^2 = F² - 2dF + d² - JuMP.add_to_expression!(ex, 1.0, F, F) - JuMP.add_to_expression!(ex, -2.0 * demand, F) - JuMP.add_to_expression!(ex, demand^2) - else - # Error term d*(1 - F/d)^2 = F²/d - 2F + d - JuMP.add_to_expression!(ex, 1.0 / demand, F, F) - JuMP.add_to_expression!(ex, -2.0, F) - JuMP.add_to_expression!(ex, demand) + # Objective function ∑ |F - d| + JuMP.set_normalized_rhs(constraint_abs_positive, -demand) + JuMP.set_normalized_rhs(constraint_abs_negative, demand) + return nothing +end + +""" +Add a term to the expression of the objective function corresponding to +the demand of a UserDemand. +""" +function add_user_demand_term!( + edge::Tuple{NodeID, NodeID}, + demand::Float64, + problem::JuMP.Model, +)::Nothing + node_id_user_demand = edge[2] + + constraint_abs_positive = problem[:abs_positive_user_demand][node_id_user_demand] + constraint_abs_negative = problem[:abs_negative_user_demand][node_id_user_demand] + + add_objective_term!(demand, constraint_abs_positive, constraint_abs_negative) +end + +""" +Add a term to the expression of the objective function corresponding to +the demand of a node with a a flow demand. +""" +function add_flow_demand_term!( + edge::Tuple{NodeID, NodeID}, + demand::Float64, + problem::JuMP.Model, +)::Nothing + node_id_flow_demand = edge[2] + + constraint_abs_positive = problem[:abs_positive_flow_demand][node_id_flow_demand] + constraint_abs_negative = problem[:abs_negative_flow_demand][node_id_flow_demand] + + add_objective_term!(demand, constraint_abs_positive, constraint_abs_negative) +end + +""" +Add a term to the expression of the objective function corresponding to +the demand of a basin. +""" +function add_basin_term!(problem::JuMP.Model, demand::Float64, node_id::NodeID)::Nothing + constraint_abs_positive = get(problem[:abs_positive_basin], node_id) + constraint_abs_negative = get(problem[:abs_negative_basin], node_id) + + if isnothing(constraint_abs_positive) + return end + + add_objective_term!(demand, constraint_abs_positive, constraint_abs_negative) return nothing end @@ -40,17 +81,29 @@ function set_objective_priority!( (; node_id, demand_reduced) = user_demand (; main_network_connections, subnetwork_demands) = allocation edge_ids = graph[].edge_ids[allocation_network_id] - F = problem[:F] - ex = JuMP.QuadExpr() + ex = JuMP.AffExpr() + + F_abs_user_demand = problem[:F_abs_user_demand] + F_abs_level_demand = problem[:F_abs_level_demand] + F_abs_flow_demand = problem[:F_abs_flow_demand] + + if !isempty(only(F_abs_user_demand.axes)) + ex += sum(F_abs_user_demand) + end + if !isempty(only(F_abs_level_demand.axes)) + ex += sum(F_abs_level_demand) + end + if !isempty(only(F_abs_flow_demand.axes)) + ex += sum(F_abs_flow_demand) + end # Terms for subnetworks as UserDemand if is_main_network(allocation_network_id) - for connections_subnetwork in main_network_connections[2:end] + for connections_subnetwork in main_network_connections for connection in connections_subnetwork d = subnetwork_demands[connection][priority_idx] - F_inlet = F[connection] - add_objective_term!(ex, d, F_inlet) + add_user_demand_term!(connection, d, problem) end end end @@ -63,8 +116,7 @@ function set_objective_priority!( # UserDemand user_demand_idx = findsorted(node_id, to_node_id) d = demand_reduced[user_demand_idx, priority_idx] - F_ud = F[edge_id] - add_objective_term!(ex, d, F_ud) + add_user_demand_term!(edge_id, d, problem) else has_demand, demand_node_id = has_external_demand(graph, to_node_id, :flow_demand) @@ -75,9 +127,8 @@ function set_objective_priority!( priority_idx == flow_priority_idx ? flow_demand.demand[findsorted(flow_demand.node_id, demand_node_id)] : 0.0 - F_fd = F[edge_id] - add_objective_term!(ex, d, F_fd) + add_flow_demand_term!(edge_id, d, problem) end end end @@ -91,8 +142,7 @@ function set_objective_priority!( get_basin_demand(allocation_model, u, p, t, node_id) : 0.0 _, basin_idx = id_index(basin.node_id, node_id) basin.demand[basin_idx] = d - F_ld = F_basin_in[node_id] - add_objective_term!(ex, d, F_ld) + add_basin_term!(problem, d, node_id) end new_objective = JuMP.@expression(problem, ex) @@ -209,8 +259,7 @@ function set_initial_capacities_source!( for edge_id in edge_ids if graph[edge_id...].allocation_network_id_source == allocation_network_id # If it is a source edge for this allocation problem - if (edge_id ∉ main_network_source_edges) | - is_main_network(allocation_network_id) + if edge_id ∉ main_network_source_edges # Reset the source to the current flow from the physical layer. source_capacity = get_flow(graph, edge_id..., 0) JuMP.set_normalized_rhs( diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 262c4e736..d804d8020 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -44,7 +44,7 @@ @test allocated[3, :] ≈ [0.0, 2.0] end -@testitem "Allocation objective" begin +@testitem "Allocation objective: linear absolute" begin using DataFrames: DataFrame using SciMLBase: successful_retcode using Ribasim: NodeID @@ -57,23 +57,15 @@ end config = Ribasim.Config(toml_path) model = Ribasim.run(config) @test successful_retcode(model) - (; p) = model.integrator - (; user_demand) = p - problem = p.allocation.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 isa JuMP.AffExpr # Affine expression + @test :F_abs_user_demand in keys(problem.obj_dict) F = problem[:F] + F_abs_user_demand = problem[:F_abs_user_demand] - to_user_5 = F[(NodeID(:Basin, 4), NodeID(:UserDemand, 5))] - to_user_6 = F[(NodeID(:Basin, 4), NodeID(:UserDemand, 6))] - - @test objective.aff.constant ≈ sum(user_demand.demand) - @test objective.aff.terms[to_user_5] ≈ -2.0 - @test objective.aff.terms[to_user_6] ≈ -2.0 - @test objective.terms[JuMP.UnorderedPair(to_user_5, to_user_5)] ≈ - 1 / user_demand.demand[1] - @test objective.terms[JuMP.UnorderedPair(to_user_6, to_user_6)] ≈ - 1 / user_demand.demand[2] + @test objective.terms[F_abs_user_demand[NodeID(:UserDemand, 5)]] == 1.0 + @test objective.terms[F_abs_user_demand[NodeID(:UserDemand, 6)]] == 1.0 end @testitem "Allocation with controlled fractional flow" begin @@ -168,7 +160,7 @@ end @test Ribasim.is_main_network(first(allocation_network_ids)) # Connections from main network to subnetworks - @test only(main_network_connections[1]) == (NodeID(:FlowBoundary, 1), NodeID(:Basin, 2)) + @test isempty(main_network_connections[1]) @test only(main_network_connections[2]) == (NodeID(:Basin, 2), NodeID(:Pump, 11)) @test only(main_network_connections[3]) == (NodeID(:Basin, 6), NodeID(:Pump, 24)) @test only(main_network_connections[4]) == (NodeID(:Basin, 10), NodeID(:Pump, 38)) @@ -181,6 +173,14 @@ end (NodeID(:Basin, 10), NodeID(:Pump, 38)), ] ⊆ allocation_edges_main_network + # Subnetworks interpreted as user_demands require variables and constraints to + # support absolute value expressions in the objective function + allocation_model_main_network = Ribasim.get_allocation_model(p, Int32(1)) + problem = allocation_model_main_network.problem + @test problem[:F_abs_user_demand].axes[1] == NodeID.(:Pump, [11, 24, 38]) + @test problem[:abs_positive_user_demand].axes[1] == NodeID.(:Pump, [11, 24, 38]) + @test problem[:abs_negative_user_demand].axes[1] == NodeID.(:Pump, [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, Int32(3)).problem[:source].axes[1] == @@ -223,40 +223,33 @@ end # See the difference between these values here and in # "subnetworks_with_sources" - @test subnetwork_demands[(NodeID(:Basin, 2), NodeID(:Pump, 11))] ≈ [4.0, 4.0, 0.0] atol = - 1e-4 + @test subnetwork_demands[(NodeID(:Basin, 2), NodeID(:Pump, 11))] ≈ [4.0, 4.0, 0.0] @test subnetwork_demands[(NodeID(:Basin, 6), NodeID(:Pump, 24))] ≈ [0.004, 0.0, 0.0] - @test subnetwork_demands[(NodeID(:Basin, 10), NodeID(:Pump, 38))][1:2] ≈ [0.001, 0.002] rtol = - 1e-4 + @test subnetwork_demands[(NodeID(:Basin, 10), NodeID(:Pump, 38))][1:2] ≈ [0.001, 0.002] # Solving for the main network, containing subnetworks as UserDemands allocation_model = allocation_models[1] (; problem) = allocation_model - Ribasim.allocate_priority!(allocation_model, u, p, t, 1, OptimizationType.allocate) + Ribasim.allocate!(p, allocation_model, t, u, OptimizationType.allocate) # Main network objective function - F = problem[:F] objective = JuMP.objective_function(problem) - objective_edges = keys(objective.terms) - F_1 = F[(NodeID(:Basin, 2), NodeID(:Pump, 11))] - F_2 = F[(NodeID(:Basin, 6), NodeID(:Pump, 24))] - F_3 = F[(NodeID(:Basin, 10), NodeID(:Pump, 38))] - @test JuMP.UnorderedPair(F_1, F_1) ∈ objective_edges - @test JuMP.UnorderedPair(F_2, F_2) ∈ objective_edges - @test JuMP.UnorderedPair(F_3, F_3) ∈ objective_edges + objective_variables = keys(objective.terms) + F_abs_user_demand = problem[:F_abs_user_demand] + @test F_abs_user_demand[NodeID(:Pump, 11)] ∈ objective_variables + @test F_abs_user_demand[NodeID(:Pump, 24)] ∈ objective_variables + @test F_abs_user_demand[NodeID(:Pump, 38)] ∈ objective_variables # Running full allocation algorithm Ribasim.set_flow!(graph, NodeID(:FlowBoundary, 1), NodeID(:Basin, 2), 4.5) u = ComponentVector(; storage = zeros(length(p.basin.node_id))) Ribasim.update_allocation!((; p, t, u)) - @test subnetwork_allocateds[NodeID(:Basin, 2), NodeID(:Pump, 11)] ≈ [4.0, 0.4947, 0.0] atol = - 1e-4 - @test subnetwork_allocateds[NodeID(:Basin, 6), NodeID(:Pump, 24)] ≈ [0.004, 0.0, 0.0] rtol = - 1e-3 - - @test subnetwork_allocateds[NodeID(:Basin, 10), NodeID(:Pump, 38)] ≈ - [0.001, 2.473e-4, 0.0] rtol = 1e-3 + @test subnetwork_allocateds[NodeID(:Basin, 2), NodeID(:Pump, 11)] ≈ + [4.0, 0.49500000, 0.0] + @test subnetwork_allocateds[NodeID(:Basin, 6), NodeID(:Pump, 24)] ≈ + [0.00399999999, 0.0, 0.0] + @test subnetwork_allocateds[NodeID(:Basin, 10), NodeID(:Pump, 38)] ≈ [0.001, 0.0, 0.0] # Test for existence of edges in allocation flow record allocation_flow = DataFrame(record_flow) @@ -269,7 +262,7 @@ end @test all(allocation_flow.edge_exists) @test user_demand.allocated[2, :] ≈ [4.0, 0.0, 0.0] - @test all(isapprox.(user_demand.allocated[7, :], [0.001, 0.0, 0.0], atol = 1e-5)) + @test user_demand.allocated[7, :] ≈ [0.001, 0.0, 0.0] end @testitem "subnetworks with sources" begin @@ -299,7 +292,7 @@ end # Collecting demands u = ComponentVector(; storage = zeros(length(basin.node_id))) - for (i, allocation_model) in enumerate(allocation_models[2:end]) + for allocation_model in allocation_models[2:end] Ribasim.allocate!(p, allocation_model, t, u, OptimizationType.internal_sources) Ribasim.allocate!(p, allocation_model, t, u, OptimizationType.collect_demands) end @@ -307,8 +300,7 @@ end # See the difference between these values here and in # "allocation with main network optimization problem", internal sources # lower the subnetwork demands - @test subnetwork_demands[(NodeID(:Basin, 2), NodeID(:Pump, 11))] ≈ [3.1, 4.0, 0.0] atol = - 1e-4 + @test subnetwork_demands[(NodeID(:Basin, 2), NodeID(:Pump, 11))] ≈ [3.1, 4.0, 0.0] @test subnetwork_demands[(NodeID(:Basin, 6), NodeID(:Pump, 24))] ≈ [0.004, 0.0, 0.0] @test subnetwork_demands[(NodeID(:Basin, 10), NodeID(:Pump, 38))][1:2] ≈ [0.001, 0.001] end @@ -407,6 +399,7 @@ end F = problem[:F] F_flow_buffer_in = problem[:F_flow_buffer_in] F_flow_buffer_out = problem[:F_flow_buffer_out] + F_abs_flow_demand = problem[:F_abs_flow_demand] node_id_with_flow_demand = NodeID(NodeType.TabulatedRatingCurve, 2) constraint_flow_out = problem[:flow_demand_outflow][node_id_with_flow_demand] @@ -442,8 +435,9 @@ end optimization_type, ) objective = JuMP.objective_function(problem) + @test F_abs_flow_demand[node_id_with_flow_demand] in keys(objective.terms) # Reduced demand - @test flow_demand.demand[1] ≈ flow_demand.demand_itp[1](t) - 0.001 rtol = 1e-3 + @test flow_demand.demand[1] == flow_demand.demand_itp[1](t) - 0.001 @test JuMP.normalized_rhs(constraint_flow_out) == Inf ## Priority 2 @@ -456,9 +450,9 @@ end optimization_type, ) # No demand left - @test flow_demand.demand[1] < 1e-10 + @test flow_demand.demand[1] ≈ 0.0 # Allocated - @test JuMP.value(only(F_flow_buffer_in)) ≈ 0.001 rtol = 1e-3 + @test JuMP.value(only(F_flow_buffer_in)) == 0.001 @test JuMP.normalized_rhs(constraint_flow_out) == 0.0 ## Priority 3 @@ -474,10 +468,9 @@ end # The flow from the source is used up in previous priorities @test JuMP.value(F[(NodeID(NodeType.LevelBoundary, 1), node_id_with_flow_demand)]) == 0 # So flow from the flow buffer is used for UserDemand #4 - @test JuMP.value(F_flow_buffer_out[node_id_with_flow_demand]) ≈ 0.001 rtol = 1e-3 - + @test JuMP.value(F_flow_buffer_out[node_id_with_flow_demand]) == 0.001 # Flow taken from buffer - @test JuMP.value(only(F_flow_buffer_out)) ≈ user_demand.demand_itp[1][3](t) rtol = 1e-3 + @test JuMP.value(only(F_flow_buffer_out)) == user_demand.demand_itp[1][3](t) # No flow coming from level boundary @test JuMP.value(F[(only(level_boundary.node_id), node_id_with_flow_demand)]) == 0 @@ -492,9 +485,8 @@ end ) # Get demand from buffers d = user_demand.demand_itp[3][4](t) - @test JuMP.value(F[(NodeID(NodeType.UserDemand, 4), NodeID(NodeType.Basin, 7))]) + - JuMP.value(F[(NodeID(NodeType.UserDemand, 6), NodeID(NodeType.Basin, 7))]) ≈ d rtol = - 1e-3 + @assert JuMP.value(F[(NodeID(NodeType.UserDemand, 4), NodeID(NodeType.Basin, 7))]) + + JuMP.value(F[(NodeID(NodeType.UserDemand, 6), NodeID(NodeType.Basin, 7))]) == d end @testitem "flow_demand_with_max_flow_rate" begin @@ -518,28 +510,3 @@ end ) @test constraint.set.upper == 2.0 end - -@testitem "equal_fraction_allocation" begin - using Ribasim: NodeID, NodeType - using StructArrays: StructVector - using DataFrames: DataFrame - - toml_path = - normpath(@__DIR__, "../../generated_testmodels/fair_distribution/ribasim.toml") - @test ispath(toml_path) - model = Ribasim.run(toml_path) - (; user_demand, graph) = model.integrator.p - - data_allocation = DataFrame(Ribasim.allocation_table(model)) - fractions = Vector{Float64}[] - - for id in user_demand.node_id - data_allocation_id = filter(:node_id => ==(id.value), data_allocation) - frac = data_allocation_id.allocated ./ data_allocation_id.demand - push!(fractions, frac) - end - - @test all(isapprox.(fractions[1], fractions[2], atol = 1e-4)) - @test all(isapprox.(fractions[1], fractions[3], atol = 1e-4)) - @test all(isapprox.(fractions[1], fractions[4], atol = 1e-4)) -end diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index b0aa5e336..e3669da7d 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -158,23 +158,25 @@ $$ The error between the flows and user demands is denoted by $E_{\text{user demand}}$, where $$ - E_{\text{user demand}} = \sum_{(i,j)\in E_S\;:\; i\in U_S} d_j^p(t)\left(1 - \frac{F_{ij}}{d_j^p(t)}\right)^2 + E_{\text{user demand}} = \sum_{(i,j)\in E_S\;:\; i\in U_S} \left| F_{ij} - d_j^p(t)\right| $$ :::{.callout-note} When performing main network allocation, the connections to subnetworks are also interpreted as UserDemand with demands determined by subnetwork demand collection. ::: -This type of objective cares about the fraction of the demand allocated, and will lead to an equal fraction of all demands allocated when possible. +This type of objective cares about the absolute amount of water allocated to a demand. It treats all deviations equally which means it doesn't give larger punishment per flow unit if deviations increase. + +The absolute value applied here is not supported in a linear programming context directly; this requires introduction of new variables and constraints. For more details see [here](https://optimization.cbe.cornell.edu/index.php?title=Optimization_with_absolute_values). Likewise, the error of level demands from basins is the absolute difference between flows consumed by basins and basin demands. $$ - E_{\text{level demand}} = \sum_{i \in B_S} d_i^p(t)\left(1 - \frac{F_i^\text{basin in}}{d_i^p(t)}\right)^2 + E_{\text{level demand}} = \sum_{i \in B_S} \left| F_i^\text{basin in} - d_i^p(t)\right| $$ Lastly, the error of the flow demands is given as below. $$ - E_{\text{flow demand}} = \sum_{i \in FD_S} d_i^p(t)\left(1 - \frac{F_i^\text{buffer in}}{d_i^p(t)}\right)^2 + E_{\text{flow demand}} = \sum_{i \in FD_S} \left| F_i^\text{buffer in} - d_i^p(t)\right| $$ ## The optimization constraints @@ -236,7 +238,7 @@ $$ {#eq-fractionalflowconstraint} ## Example -The following is an example of an optimization problem for the example model shown [here](../python/examples.ipynb#model-with-allocation-user-demand): +The following is an example of an optimization problem for the example shown [here](../python/examples.ipynb#model-with-allocation-user-demand): ```{julia} diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index 3706e2173..ceeb1fad2 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -7,7 +7,6 @@ import ribasim_testmodels from ribasim_testmodels.allocation import ( allocation_example_model, - fair_distribution_model, flow_demand_model, fractional_flow_subnetwork_model, level_demand_model, @@ -67,7 +66,6 @@ "bucket_model", "discrete_control_of_pid_control_model", "dutch_waterways_model", - "fair_distribution_model", "flow_boundary_time_model", "flow_condition_model", "flow_demand_model", diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 5172a8dd8..1e102e8c0 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -1029,100 +1029,3 @@ def linear_resistance_demand_model(): model.edge.add(model.flow_demand[4], model.linear_resistance[2]) return model - - -def fair_distribution_model(): - model = Model( - starttime="2020-01-01 00:00:00", - endtime="2020-01-07 00:00:00", - crs="EPSG:28992", - allocation=Allocation(use_allocation=True), - ) - - model.level_boundary.add( - Node(1, Point(0, 0), subnetwork_id=1), - [ - level_boundary.Static( - level=[1.0], - ) - ], - ) - - model.pump.add( - Node( - 2, - Point(1, 0), - subnetwork_id=1, - ), - [pump.Static(flow_rate=9.0, max_flow_rate=[9.0])], - ) - - model.basin.add( - Node(3, Point(2, 0), subnetwork_id=1), - [basin.Profile(area=1e3, level=[0.0, 1.0]), basin.State(level=[1.0])], - ) - - model.linear_resistance.add( - Node(4, Point(3, 0), subnetwork_id=1), - [linear_resistance.Static(resistance=[1.0])], - ) - - model.basin.add( - Node(5, Point(4, 0), subnetwork_id=1), - [basin.Profile(area=1e3, level=[0.0, 1.0]), basin.State(level=[1.0])], - ) - - model.user_demand.add( - Node(6, Point(2, 1), subnetwork_id=1), - [ - user_demand.Static( - priority=[1], demand=1.0, return_factor=1.0, min_level=0.2 - ) - ], - ) - - model.user_demand.add( - Node(7, Point(2, -1), subnetwork_id=1), - [ - user_demand.Static( - priority=[1], demand=2.0, return_factor=1.0, min_level=0.2 - ) - ], - ) - - model.user_demand.add( - Node(8, Point(4, 1), subnetwork_id=1), - [ - user_demand.Static( - priority=[1], demand=3.0, return_factor=1.0, min_level=0.2 - ) - ], - ) - - model.user_demand.add( - Node(9, Point(4, -1), subnetwork_id=1), - [ - user_demand.Time( - priority=1, - time=pd.date_range(start="2020-01", end="2021-01", freq="MS"), - demand=np.linspace(1.0, 5.0, 13), - return_factor=1.0, - min_level=0.2, - ) - ], - ) - - model.edge.add(model.level_boundary[1], model.pump[2], subnetwork_id=1) - model.edge.add(model.pump[2], model.basin[3]) - model.edge.add(model.basin[3], model.linear_resistance[4]) - model.edge.add(model.linear_resistance[4], model.basin[5]) - model.edge.add(model.basin[3], model.user_demand[6]) - model.edge.add(model.basin[3], model.user_demand[7]) - model.edge.add(model.basin[5], model.user_demand[8]) - model.edge.add(model.basin[5], model.user_demand[9]) - model.edge.add(model.user_demand[6], model.basin[3]) - model.edge.add(model.user_demand[7], model.basin[3]) - model.edge.add(model.user_demand[8], model.basin[5]) - model.edge.add(model.user_demand[9], model.basin[5]) - - return model