Revert "Equal ratio allocation objective (#1366)"

This reverts commit 19139f1.

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,13 +921,18 @@ 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,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(