Skip to content

Commit

Permalink
Equal ratio allocation objective (#1366)
Browse files Browse the repository at this point in the history 
Fixes #1348.

The objective function I found that leads to equal fraction allocation
has terms of the form

$$
d\left(1 - \frac{F}{d}\right)^2 = \left(\sqrt{d} -
\frac{F}{\sqrt{d}}\right)^2 = \frac{F^2}{d} - 2F + d.
$$

I can't explain yet why this works.
  • Loading branch information
@SouthEndMusic
SouthEndMusic authored Apr 12, 2024
1 parent 73ab094 commit 19139f1
Show file tree
Hide file tree
Showing 6 changed files with 208 additions and 281 deletions.
158 changes: 2 additions & 156 deletions core/src/allocation_init.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,8 @@ 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
if (graph[outflow_id].allocation_network_id != 1) |
(graph[node_id, outflow_id].allocation_network_id_source == 1)
main_network_source_edges =
get_main_network_connections(p, graph[outflow_id].allocation_network_id)
edge = (node_id, outflow_id)
Expand Down Expand Up @@ -405,55 +406,6 @@ 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.
Expand All @@ -468,7 +420,6 @@ 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}[]
Expand Down Expand Up @@ -692,106 +643,6 @@ 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.
Expand Down Expand Up @@ -921,18 +772,13 @@ 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
add_constraints_conservation_basin!(problem, p, allocation_network_id)
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)
Expand Down
105 changes: 28 additions & 77 deletions core/src/allocation_optim.jl
View file
Original file line number Diff line number Diff line change
@@ -1,67 +1,26 @@
@enumx OptimizationType internal_sources collect_demands allocate

"""
Add a term to the objective function given by the objective type,
depending in the provided flow variable and the associated demand.
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.
"""
function add_objective_term!(
ex::JuMP.QuadExpr,
demand::Float64,
constraint_abs_positive::Union{JuMP.ConstraintRef, Nothing} = nothing,
constraint_abs_negative::Union{JuMP.ConstraintRef, Nothing} = nothing,
F::JuMP.VariableRef,
)::Nothing
# 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
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)
end

add_objective_term!(demand, constraint_abs_positive, constraint_abs_negative)
return nothing
end

Expand All @@ -81,29 +40,17 @@ 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.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
ex = JuMP.QuadExpr()

# Terms for subnetworks as UserDemand
if is_main_network(allocation_network_id)
for connections_subnetwork in main_network_connections
for connections_subnetwork in main_network_connections[2:end]
for connection in connections_subnetwork
d = subnetwork_demands[connection][priority_idx]
add_user_demand_term!(connection, d, problem)
F_inlet = F[connection]
add_objective_term!(ex, d, F_inlet)
end
end
end
Expand All @@ -116,7 +63,8 @@ function set_objective_priority!(
# UserDemand
user_demand_idx = findsorted(node_id, to_node_id)
d = demand_reduced[user_demand_idx, priority_idx]
add_user_demand_term!(edge_id, d, problem)
F_ud = F[edge_id]
add_objective_term!(ex, d, F_ud)
else
has_demand, demand_node_id =
has_external_demand(graph, to_node_id, :flow_demand)
Expand All @@ -127,8 +75,9 @@ 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_flow_demand_term!(edge_id, d, problem)
add_objective_term!(ex, d, F_fd)
end
end
end
Expand All @@ -142,7 +91,8 @@ 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!(problem, d, node_id)
F_ld = F_basin_in[node_id]
add_objective_term!(ex, d, F_ld)
end

new_objective = JuMP.@expression(problem, ex)
Expand Down Expand Up @@ -259,7 +209,8 @@ 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
if (edge_id ∉ main_network_source_edges) |
is_main_network(allocation_network_id)
# Reset the source to the current flow from the physical layer.
source_capacity = get_flow(graph, edge_id..., 0)
JuMP.set_normalized_rhs(
Expand Down
Loading

0 comments on commit 19139f1

Please sign in to comment.