diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 259a30fea..795e83d4f 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -399,38 +399,36 @@ function add_variables_absolute_value!( problem::JuMP.Model, p::Parameters, allocation_network_id::Int, - config::Config, )::Nothing (; 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_demand = NodeID[] - node_ids_basin = 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_basin, node_id) - end + + node_ids = graph[].node_ids[allocation_network_id] + node_ids_user_demand = NodeID[] + node_ids_basin = 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_basin, 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 + # 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 - - problem[:F_abs_user_demand] = - JuMP.@variable(problem, F_abs_user_demand[node_id = node_ids_user_demand]) - problem[:F_abs_basin] = - JuMP.@variable(problem, F_abs_basin[node_id = node_ids_basin]) end + + problem[:F_abs_user_demand] = + JuMP.@variable(problem, F_abs_user_demand[node_id = node_ids_user_demand]) + problem[:F_abs_basin] = JuMP.@variable(problem, F_abs_basin[node_id = node_ids_basin]) + return nothing end @@ -632,7 +630,6 @@ function add_constraints_absolute_value!( problem::JuMP.Model, flow_per_node::Dict{NodeID, JuMP.VariableRef}, F_abs::JuMP.Containers.DenseAxisArray, - objective_type::String, variable_type::String, )::Nothing # Example demand @@ -640,41 +637,23 @@ function add_constraints_absolute_value!( node_ids = only(F_abs.axes) - if objective_type == "linear_absolute" - # 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 - ) - elseif objective_type == "linear_relative" - # These constraints together make sure that F_abs_user_demand acts as the absolute - # value F_abs_user_demand = |x| where x = 1-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] >= (1 - 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] >= -(1 - flow_per_node[node_id] / d), - base_name = base_name - ) - end + # 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 @@ -685,28 +664,24 @@ 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, - config::Config, )::Nothing (; graph) = p - objective_type = config.allocation.objective_type - if startswith(objective_type, "linear") - F = problem[:F] - F_abs_user_demand = problem[:F_abs_user_demand] + 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) - ) + 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", + ) - add_constraints_absolute_value!( - problem, - flow_per_node, - F_abs_user_demand, - objective_type, - "user_demand", - ) - end return nothing end @@ -714,22 +689,14 @@ end Add constraints so that variables F_abs_basin act as the absolute value of the expression comparing flow to a basin to its demand. """ -function add_constraints_absolute_value_basin!(problem::JuMP.Model, config::Config)::Nothing - objective_type = config.allocation.objective_type - if startswith(objective_type, "linear") - F_basin_in = problem[:F_basin_in] - F_abs_basin = problem[:F_abs_basin] - flow_per_node = - Dict(node_id => F_basin_in[node_id] for node_id in only(F_abs_basin.axes)) +function add_constraints_absolute_value_basin!(problem::JuMP.Model)::Nothing + F_basin_in = problem[:F_basin_in] + F_abs_basin = problem[:F_abs_basin] + flow_per_node = + Dict(node_id => F_basin_in[node_id] for node_id in only(F_abs_basin.axes)) + + add_constraints_absolute_value!(problem, flow_per_node, F_abs_basin, "basin") - add_constraints_absolute_value!( - problem, - flow_per_node, - F_abs_basin, - objective_type, - "basin", - ) - end return nothing end @@ -806,7 +773,6 @@ end Construct the allocation problem for the current subnetwork as a JuMP.jl model. """ function allocation_problem( - config::Config, p::Parameters, capacity::SparseMatrixCSC{Float64, Int}, allocation_network_id::Int, @@ -817,15 +783,15 @@ 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, config) + add_variables_absolute_value!(problem, p, allocation_network_id) # Add constraints to problem add_constraints_capacity!(problem, capacity, p, allocation_network_id) add_constraints_source!(problem, p, allocation_network_id) add_constraints_flow_conservation!(problem, p, allocation_network_id) add_constraints_user_demand_returnflow!(problem, p, allocation_network_id) - add_constraints_absolute_value_user_demand!(problem, p, config) - add_constraints_absolute_value_basin!(problem, config) + add_constraints_absolute_value_user_demand!(problem, p) + add_constraints_absolute_value_basin!(problem) add_constraints_fractional_flow!(problem, p, allocation_network_id) add_constraints_basin_flow!(problem) @@ -837,7 +803,6 @@ Construct the JuMP.jl problem for allocation. Inputs ------ -config: The model configuration with allocation configuration in config.allocation p: Ribasim problem parameters Δt_allocation: The timestep between successive allocation solves @@ -846,7 +811,6 @@ Outputs An AllocationModel object. """ function AllocationModel( - config::Config, allocation_network_id::Int, p::Parameters, Δt_allocation::Float64, @@ -855,16 +819,7 @@ function AllocationModel( capacity = allocation_graph(p, allocation_network_id) # The JuMP.jl allocation problem - problem = allocation_problem(config, p, capacity, allocation_network_id) - if config.allocation.objective_type != "linear_absolute" - error("Type of object function is not supported") - end + problem = allocation_problem(p, capacity, allocation_network_id) - return AllocationModel( - Symbol(config.allocation.objective_type), - allocation_network_id, - capacity, - problem, - Δt_allocation, - ) + return AllocationModel(allocation_network_id, capacity, problem, Δt_allocation) end diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 9593b68e4..ac5c00d08 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -3,48 +3,13 @@ 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::Union{JuMP.QuadExpr, JuMP.AffExpr}, - F_variable::JuMP.VariableRef, demand::Float64, - objective_type::Symbol; constraint_abs_positive::Union{JuMP.ConstraintRef, Nothing} = nothing, constraint_abs_negative::Union{JuMP.ConstraintRef, Nothing} = nothing, )::Nothing - if objective_type == :quadratic_absolute - # Objective function ∑ (F - d)^2 - JuMP.add_to_expression!(ex, 1, F_variable, F_variable) - JuMP.add_to_expression!(ex, -2 * demand, F_variable) - 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_variable, F_variable) - JuMP.add_to_expression!(ex, -2.0 / demand, F_variable) - JuMP.add_to_expression!(ex, 1.0) - - elseif objective_type == :linear_absolute - # Objective function ∑ |F - d| - JuMP.set_normalized_rhs(constraint_abs_positive, -demand) - JuMP.set_normalized_rhs(constraint_abs_negative, demand) - - elseif objective_type == :linear_relative - # Objective function ∑ |1 - F/d| - JuMP.set_normalized_coefficient( - constraint_abs_positive, - F_variable, - iszero(demand) ? 0 : 1 / demand, - ) - JuMP.set_normalized_coefficient( - constraint_abs_negative, - F_variable, - iszero(demand) ? 0 : -1 / demand, - ) - else - error("Invalid allocation objective type $objective_type.") - end + # Objective function ∑ |F - d| + JuMP.set_normalized_rhs(constraint_abs_positive, -demand) + JuMP.set_normalized_rhs(constraint_abs_negative, demand) return nothing end @@ -53,64 +18,27 @@ Add a term to the expression of the objective function corresponding to the demand of a UserDemand. """ function add_user_demand_term!( - ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, edge::Tuple{NodeID, NodeID}, - objective_type::Symbol, demand::Float64, problem::JuMP.Model, )::Nothing - F = problem[:F] - F_edge = F[edge] node_id_user_demand = edge[2] - if objective_type in [:linear_absolute, :linear_relative] - constraint_abs_positive = problem[:abs_positive_user_demand][node_id_user_demand] - constraint_abs_negative = problem[:abs_negative_user_demand][node_id_user_demand] - else - constraint_abs_positive = nothing - constraint_abs_negative = nothing - end + 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!( - ex, - F_edge, - demand, - objective_type; - constraint_abs_positive, - constraint_abs_negative, - ) + 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!( - ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, - problem::JuMP.Model, - demand::Float64, - objective_type::Symbol, - node_id::NodeID, -)::Nothing - F_basin_in = problem[:F_basin_in] - F_basin = F_basin_in[node_id] +function add_basin_term!(problem::JuMP.Model, demand::Float64, node_id::NodeID)::Nothing + constraint_abs_positive = problem[:abs_positive_basin][node_id] + constraint_abs_negative = problem[:abs_negative_basin][node_id] - if objective_type in [:linear_absolute, :linear_relative] - constraint_abs_positive = problem[:abs_positive_basin][node_id] - constraint_abs_negative = problem[:abs_negative_basin][node_id] - else - constraint_abs_positive = nothing - constraint_abs_negative = nothing - end - - add_objective_term!( - ex, - F_basin, - demand, - objective_type; - constraint_abs_positive, - constraint_abs_negative, - ) + add_objective_term!(demand, constraint_abs_positive, constraint_abs_negative) return nothing end @@ -125,19 +53,15 @@ function set_objective_priority!( t::Float64, priority_idx::Int, )::Nothing - (; objective_type, problem, allocation_network_id) = allocation_model + (; problem, allocation_network_id) = allocation_model (; graph, user_demand, allocation, basin) = p (; demand_itp, demand_from_timeseries, node_id) = user_demand (; main_network_connections, subnetwork_demands) = allocation edge_ids = graph[].edge_ids[allocation_network_id] - if objective_type in [:quadratic_absolute, :quadratic_relative] - ex = JuMP.QuadExpr() - elseif objective_type in [:linear_absolute, :linear_relative] - ex = JuMP.AffExpr() - ex += sum(problem[:F_abs_user_demand]) - ex += sum(problem[:F_abs_basin]) - end + ex = JuMP.AffExpr() + ex += sum(problem[:F_abs_user_demand]) + ex += sum(problem[:F_abs_basin]) demand_max = 0.0 @@ -147,7 +71,7 @@ function set_objective_priority!( for connection in connections_subnetwork d = subnetwork_demands[connection][priority_idx] demand_max = max(demand_max, d) - add_user_demand_term!(ex, connection, objective_type, d, problem) + add_user_demand_term!(connection, d, problem) end end end @@ -167,7 +91,7 @@ function set_objective_priority!( d = get_user_demand(p, node_id_user_demand, priority_idx) end demand_max = max(demand_max, d) - add_user_demand_term!(ex, edge_id, objective_type, d, problem) + add_user_demand_term!(edge_id, d, problem) end # Terms for basins @@ -179,7 +103,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 - add_basin_term!(ex, problem, d, objective_type, node_id) + add_basin_term!(problem, d, node_id) end new_objective = JuMP.@expression(problem, ex) diff --git a/core/src/config.jl b/core/src/config.jl index 56fbadc34..1c0020369 100644 --- a/core/src/config.jl +++ b/core/src/config.jl @@ -115,7 +115,6 @@ end @option struct Allocation <: TableOption timestep::Union{Float64, Nothing} = nothing use_allocation::Bool = false - objective_type::String = "linear_absolute" end @option @addnodetypes struct Toml <: TableOption diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 6d02cbceb..6e8091bf6 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -43,14 +43,12 @@ const VectorInterpolation = """ Store information for a subnetwork used for allocation. -objective_type: The name of the type of objective used allocation_network_id: The ID of this allocation network capacity: The capacity per edge of the allocation network, as constrained by nodes that have a max_flow_rate problem: The JuMP.jl model for solving the allocation problem Δt_allocation: The time interval between consecutive allocation solves """ struct AllocationModel - objective_type::Symbol allocation_network_id::Int capacity::SparseMatrixCSC{Float64, Int} problem::JuMP.Model diff --git a/core/src/read.jl b/core/src/read.jl index 4e4bdcf55..dec5df4b8 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -226,7 +226,7 @@ function initialize_allocation!(p::Parameters, config::Config)::Nothing for allocation_network_id in allocation_network_ids_ push!( allocation_models, - AllocationModel(config, allocation_network_id, p, config.allocation.timestep), + AllocationModel(allocation_network_id, p, config.allocation.timestep), ) end return nothing diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 985ab1cc8..c0e45b60d 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -49,34 +49,6 @@ @test Ribasim.get_user_demand(p, NodeID(:UserDemand, 11), 2) ≈ π end -@testitem "Allocation objective: quadratic absolute" begin - using DataFrames: DataFrame - using SciMLBase: successful_retcode - using Ribasim: NodeID - import JuMP - - toml_path = - normpath(@__DIR__, "../../generated_testmodels/minimal_subnetwork/ribasim.toml") - @test ispath(toml_path) - - config = Ribasim.Config(toml_path; allocation_objective_type = "quadratic_absolute") - @test_throws "Type of object function is not supported" model = Ribasim.run(config) -end - -@testitem "Allocation objective: quadratic relative" begin - using DataFrames: DataFrame - using SciMLBase: successful_retcode - using Ribasim: NodeID - import JuMP - - toml_path = - normpath(@__DIR__, "../../generated_testmodels/minimal_subnetwork/ribasim.toml") - @test ispath(toml_path) - - config = Ribasim.Config(toml_path; allocation_objective_type = "quadratic_relative") - @test_throws "Type of object function is not supported" model = Ribasim.run(config) -end - @testitem "Allocation objective: linear absolute" begin using DataFrames: DataFrame using SciMLBase: successful_retcode @@ -87,7 +59,7 @@ end normpath(@__DIR__, "../../generated_testmodels/minimal_subnetwork/ribasim.toml") @test ispath(toml_path) - config = Ribasim.Config(toml_path; allocation_objective_type = "linear_absolute") + config = Ribasim.Config(toml_path) model = Ribasim.run(config) @test successful_retcode(model) problem = model.integrator.p.allocation.allocation_models[1].problem @@ -101,20 +73,6 @@ end @test objective.terms[F_abs_user_demand[NodeID(:UserDemand, 6)]] == 1.0 end -@testitem "Allocation objective: linear relative" begin - using DataFrames: DataFrame - using SciMLBase: successful_retcode - using Ribasim: NodeID - import JuMP - - toml_path = - normpath(@__DIR__, "../../generated_testmodels/minimal_subnetwork/ribasim.toml") - @test ispath(toml_path) - - config = Ribasim.Config(toml_path; allocation_objective_type = "linear_relative") - @test_throws "Type of object function is not supported" model = Ribasim.run(config) -end - @testitem "Allocation with controlled fractional flow" begin using DataFrames using Ribasim: NodeID diff --git a/core/test/docs.toml b/core/test/docs.toml index d235522f6..60b98a43e 100644 --- a/core/test/docs.toml +++ b/core/test/docs.toml @@ -18,7 +18,6 @@ time = "basin/time.arrow" [allocation] timestep = 86400 # optional (required if use_allocation = true), default 86400 use_allocation = false # optional, default false -objective_type = "linear_absolute" # optional, default "linear_absolute" [solver] algorithm = "QNDF" # optional, default "QNDF" diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index 73be116e9..681295e66 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -122,69 +122,29 @@ There are several types of variable whose value has to be determined to solve th ## The optimization objective -The goal of allocation is to get the flow to nodes with demands as close as possible to these demands. To achieve this, a sum error of terms is minimized. The form of these error terms is determined by the choice of one of the supported objective function types, which are discussed further below. +The goal of allocation is to get the flow to nodes with demands as close as possible to these demands. To achieve this, a sum error of terms is minimized. $$ - \min E_{\text{user_demand}} + E_{\text{basin}} + \min E_{\text{user demand}} + E_{\text{basin}} $$ -### User demands - -- `quadratic_absolute`: -$$ - E_{\text{user_demand}} = \sum_{(i,j)\in E_S\;:\; i\in U_S} \left( F_{ij} - d_j^p(t)\right)^2 -$$ -- `quadratic_relative`: -$$ - E_{\text{user_demand}} = \sum_{(i,j)\in E_S\;:\; i\in U_S} \left( 1 - \frac{F_{ij}}{d_j^p(t)}\right)^2 -$$ -- `linear_absolute` (default): +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} \left| F_{ij} - d_j^p(t)\right| -$$ -- `linear_relative`: -$$ - E_{\text{user_demand}} = \sum_{(i,j)\in E_S\;:\; i\in U_S} \left|1 - \frac{F_{ij}}{d_j^p(t)}\right| + 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. ::: - -To avoid division by $0$ errors, if a `*_relative` objective is used and a demand is $0$, the coefficient of the flow $F_{ij}$ is set to $0$. - -For `*_absolute` objectives the optimizer cares about the actual amount of water allocated to a demand, for `*_relative` objectives it cares about the fraction of the demand allocated. -For `quadratic_*` objectives the optimizer cares about avoiding large shortages, for `linear_*` objectives it treats all deviations equally. - -:::{.callout-note} -These options for objectives for allocation to demands have not been tested thoroughly, and might change in the future. -::: +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). -:::{.callout-note} -In the future new optimization objectives will be introduced, for demands of basins and priorities over sources. These will be used in combination with the above, in the form of goal programming. -::: - -### Basin demands - -- `quadratic_absolute`: -$$ - E_{\text{basin}} = \sum_{i \in B_S} \left( F_i^\text{basin in} - d_i^p(t)\right)^2 -$$ -- `quadratic_relative`: -$$ - E_{\text{basin}} = \sum_{i \in B_S} \left( 1 - \frac{F_i^\text{basin in}}{d_i^p(t)}\right)^2 -$$ -- `linear_absolute` (default): +Likewise, the error of basin demands is the absolute difference between flows consumed by basins and basin demands. $$ E_{\text{basin}} = \sum_{i \in B_S} \left| F_i^\text{basin in} - d_i^p(t)\right| $$ -- `linear_relative`: -$$ - E_{\text{basin}} = \sum_{i \in B_S} \left|1 - \frac{F_i^\text{basin in}}{d_i^p(t)}\right| -$$ ## The optimization constraints - Flow conservation: For the basins in the allocation network we have that diff --git a/docs/python/examples.ipynb b/docs/python/examples.ipynb index 3f7573a2e..f2ab8ce82 100644 --- a/docs/python/examples.ipynb +++ b/docs/python/examples.ipynb @@ -1844,7 +1844,6 @@ "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", diff --git a/python/ribasim/ribasim/config.py b/python/ribasim/ribasim/config.py index 02f4c03c6..e4a40134e 100644 --- a/python/ribasim/ribasim/config.py +++ b/python/ribasim/ribasim/config.py @@ -38,7 +38,6 @@ class Allocation(ChildModel): timestep: float | None = None use_allocation: bool = False - objective_type: str = "linear_absolute" class Results(ChildModel):