Equal ratio allocation objective

core/src/allocation_init.jl

@@ -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)
@@ -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.
@@ -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}[]
@@ -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.
@@ -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)

core/src/allocation_optim.jl

@@ -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 (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
+        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
+        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
 
@@ -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
@@ -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)
@@ -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
@@ -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)
@@ -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(