Flow communication between primary and secondary networks

core/src/allocation_init.jl

@@ -131,17 +131,13 @@ function add_conservation!(
     flow = problem[:flow]
     inflow_sum = Dict(
         basin_id => sum(
-            flow[(other_id, basin_id)] for
-            other_id in basin.inflow_ids[basin_id.idx] if
-            graph[other_id].subnetwork_id == subnetwork_id;
+            flow[(other_id, basin_id)] for other_id in basin.inflow_ids[basin_id.idx];
             init = 0,
         ) for basin_id in basin_ids_subnetwork
     )
     outflow_sum = Dict(
         basin_id => sum(
-            flow[(basin_id, other_id)] for
-            other_id in basin.outflow_ids[basin_id.idx] if
-            graph[other_id].subnetwork_id == subnetwork_id;
+            flow[(basin_id, other_id)] for other_id in basin.outflow_ids[basin_id.idx];
             init = 0,
         ) for basin_id in basin_ids_subnetwork
     )
@@ -197,7 +193,8 @@ function add_user_demand!(
         flow[inflow_link[node_id.idx].link] == sum(
             user_demand_allocated[node_id, demand_priority] for
             demand_priority in DemandPriorityIterator(node_id, p_independent)
-        )
+        );
+        base_name = "user_demand_allocated_sum_constraint"
     )
 
     # Define decision variables: Per UserDemand node the allocation error per priority for
@@ -556,7 +553,7 @@ function add_outlet!(
     return nothing
 end
 
-function add_subnetwork_demand!(
+function add_secondary_network_demand!(
     allocation_model::AllocationModel,
     p_independent::ParametersIndependent,
 )::Nothing
@@ -568,34 +565,55 @@ function add_subnetwork_demand!(
     connecting_links =
         vcat(sort!(collect(values(allocation.primary_network_connections)))...)
 
-    # Define parameters: flow allocated to user subnetworks (scaling.flow * m^3/s, values to be filled in before optimizing)
-    subnetwork_allocated =
-        problem[:subnetwork_allocated] =
-            JuMP.@variable(problem, subnetwork[connecting_links] == 0)
-
-    # Define decision variables: lower and upper user demand error (unitless)
-    relative_subnetwork_error_lower =
-        problem[:relative_subnetwork_error_lower] =
-            JuMP.@variable(problem, relative_subnetwork_error_lower[connecting_links] ≥ 0)
-    relative_subnetwork_error_upper =
-        problem[:relative_subnetwork_error_upper] =
-            JuMP.@variable(problem, relative_subnetwork_error_upper[connecting_links] ≥ 0)
-
-    # Define constraints: error terms
+    # Define decision variables: flow allocated to secondary networks
+    # (scaling.flow * m^3/s, values to be filled in before optimizing)
     d = 2.0 # Example demand (scaling.flow * m^3/s, values to be filled in before optimizing)
-    problem[:subnetwork_constraint_lower] = JuMP.@constraint(
+    secondary_network_allocated =
+        problem[:secondary_network_allocated] = JuMP.@variable(
+            problem,
+            0 ≤
+            secondary_network_allocated[
+                link = connecting_links,
+                DemandPriorityIterator(link, p_independent),
+            ] ≤
+            d
+        )
+
+    # Define constraints: The sum of the flows allocated to the secondary network is equal to the total flow into the secondary network
+    problem[:secondary_network_allocated_sum_constraint] = JuMP.@constraint(
         problem,
         [link = connecting_links],
-        d * (relative_subnetwork_error_lower[link]) ≥
-        -(flow[link] - subnetwork_allocated[link]),
-        base_name = "subnetwork_constraint_lower"
+        flow[link] == sum(
+            secondary_network_allocated[link, demand_priority] for
+            demand_priority in DemandPriorityIterator(link, p_independent)
+        );
+        base_name = "secondary_network_allocated_sum_constraint",
     )
-    problem[:subnetwork_constraint_upper] = JuMP.@constraint(
+
+    # Define decision variables: Per secondary network connection, the allocation error per priority for
+    # which the secondary network has a demand, for both the first and second objective (unitless)
+    secondary_network_error =
+        problem[:secondary_network_error] = JuMP.@variable(
+            problem,
+            0 ≤
+            secondary_network_error[
+                link = connecting_links,
+                DemandPriorityIterator(link, p_independent),
+                [:first, :second],
+            ] ≤
+            1
+        )
+
+    # Define constraints: error terms
+    problem[:secondary_network_relative_error_constraint] = JuMP.@constraint(
         problem,
-        [link = connecting_links],
-        d * (relative_subnetwork_error_upper[link]) ≥
-        flow[link] - subnetwork_allocated[link],
-        base_name = "subnetwork_constraint_upper"
+        [
+            link = connecting_links,
+            demand_priority = DemandPriorityIterator(link, p_independent),
+        ],
+        d * secondary_network_error[link, demand_priority, :first] ≥
+        d - secondary_network_allocated[link, demand_priority];
+        base_name = "secondary_network_relative_error_constraint"
     )
 
     return nothing
@@ -609,7 +627,7 @@ function add_demand_objectives!(
     allocation_model::AllocationModel,
     p_independent::ParametersIndependent,
 )::Nothing
-    (; objectives, problem, node_ids_in_subnetwork) = allocation_model
+    (; objectives, problem, node_ids_in_subnetwork, subnetwork_id) = allocation_model
     (;
         user_demand_ids_subnetwork,
         node_ids_subnetwork_with_flow_demand,
@@ -640,19 +658,27 @@ function add_demand_objectives!(
         # Objective for a fair distribution of what was allocated with the previous objective
         second_objective_expression = JuMP.AffExpr()
 
-        for error_collection in [user_demand_error, flow_demand_error]
-            for (node_id, demand_priority_, objective_ord) in keys(error_collection.data)
+        # Add UserDemand, FlowDemand and secondary network errors
+        error_collections = [user_demand_error, flow_demand_error]
+        if is_primary_network(subnetwork_id)
+            subnetwork_demand_error = problem[:secondary_network_error]
+            push!(error_collections, subnetwork_demand_error)
+        end
+
+        for error_collection in error_collections
+            for (identifier, demand_priority_, objective_ord) in keys(error_collection.data)
                 if demand_priority == demand_priority_
                     JuMP.add_to_expression!(
                         (objective_ord == :first) ? first_objective_expression :
                         second_objective_expression,
-                        error_collection[node_id, demand_priority, objective_ord],
+                        error_collection[identifier, demand_priority, objective_ord],
                     )
                     has_flow_unit_demands = true
                 end
             end
         end
 
+        # Add LevelDemand errors
         for (node_id, demand_priority_, side, objective_ord) in
             keys(level_demand_error.data)
             if demand_priority == demand_priority_
@@ -738,9 +764,30 @@ function add_demand_objectives!(
         ],
         flow_demand_error[node_id, demand_priority, :second] ≥
         flow_demand_error[node_id, demand_priority, :first] -
-        average_flow_unit_error[demand_priority]
+        average_flow_unit_error[demand_priority],
+        base_name = "flow_demand_fairness_error_constraint"
     )
 
+    # If this is the primary network, also add fairness constraints for secondary network demands
+    if is_primary_network(subnetwork_id)
+        subnetwork_demand_error = problem[:secondary_network_error]
+        # Sort connections for deterministic problem generation
+        connecting_links =
+            vcat(sort!(collect(values(allocation.primary_network_connections)))...)
+
+        problem[:subnetwork_demand_fairness_error_constraint] = JuMP.@constraint(
+            problem,
+            [
+                link = connecting_links,
+                demand_priority = DemandPriorityIterator(link, p_independent),
+            ],
+            subnetwork_demand_error[link, demand_priority, :second] ≥
+            subnetwork_demand_error[link, demand_priority, :first] -
+            average_flow_unit_error[demand_priority],
+            base_name = "subnetwork_demand_fairness_error_constraint"
+        )
+    end
+
     # Define variables: average level error for storage unit demands (LevelDemand) per demand priority (m)
     average_storage_unit_error =
         problem[:average_storage_unit_error] = JuMP.@variable(
@@ -894,22 +941,65 @@ function NodeIDsInSubnetwork(
     return node_ids_in_subnetwork
 end
 
+function has_demand_priority_subnetwork(
+    p_independent::ParametersIndependent,
+    node_ids_in_subnetwork::NodeIDsInSubnetwork,
+)::Vector{Bool}
+    (; allocation, graph, user_demand, flow_demand, level_demand) = p_independent
+    (; demand_priorities_all) = allocation
+    (;
+        user_demand_ids_subnetwork,
+        node_ids_subnetwork_with_flow_demand,
+        level_demand_ids_subnetwork,
+    ) = node_ids_in_subnetwork
+
+    has_demand_priority = zeros(Bool, length(demand_priorities_all))
+
+    for node_id in user_demand_ids_subnetwork
+        has_demand_priority .|= view(user_demand.has_demand_priority, node_id.idx, :)
+    end
+
+    for node_id in node_ids_subnetwork_with_flow_demand
+        has_demand_priority .|= view(flow_demand.has_demand_priority, node_id.idx, :)
+    end
+
+    for node_id in level_demand_ids_subnetwork
+        has_demand_priority .|= view(level_demand.has_demand_priority, node_id.idx, :)
+    end
+
+    return has_demand_priority
+end
+
 function AllocationModel(
     subnetwork_id::Int32,
     p_independent::ParametersIndependent,
     allocation_config::config.Allocation,
 )
     Δt_allocation = allocation_config.timestep
-    problem = JuMP.Model(() -> MOA.Optimizer(get_optimizer()))
-    set_multi_objective_attributes!(problem)
+    problem = JuMP.Model()
+    JuMP.set_optimizer(problem, get_optimizer())
     node_ids_in_subnetwork = NodeIDsInSubnetwork(p_independent, subnetwork_id)
     scaling = ScalingFactors(p_independent, subnetwork_id, Δt_allocation)
+    has_demand_priority =
+        has_demand_priority_subnetwork(p_independent, node_ids_in_subnetwork)
+
+    # Initialize secondary_network_demand before constructing AllocationModel
+    secondary_network_demand = Dict{Tuple{NodeID, NodeID}, Vector{Float64}}()
+    if !is_primary_network(subnetwork_id)
+        n_priorities = length(p_independent.allocation.demand_priorities_all)
+        for link in p_independent.allocation.primary_network_connections[subnetwork_id]
+            secondary_network_demand[link] = zeros(n_priorities)
+        end
+    end
+
     allocation_model = AllocationModel(;
         subnetwork_id,
         node_ids_in_subnetwork,
         problem,
         Δt_allocation,
         scaling,
+        has_demand_priority,
+        secondary_network_demand,
     )
 
     # Volume and flow
@@ -928,21 +1018,19 @@ function AllocationModel(
     add_pump!(allocation_model, p_independent)
     add_outlet!(allocation_model, p_independent)
 
-    # Demand nodes and subnetworks as demand nodes
+    # Demand nodes
     add_user_demand!(allocation_model, p_independent)
     add_flow_demand!(allocation_model, p_independent)
     add_level_demand!(allocation_model, p_independent)
 
     # Primary to secondary subnetwork connections
     if is_primary_network(subnetwork_id)
-        add_subnetwork_demand!(allocation_model, p_independent)
-    else
-        # Initialize subnetwork demands
-        n_demands = length(p_independent.allocation.demand_priorities_all)
-        if !is_primary_network(subnetwork_id)
-            for link in p_independent.allocation.primary_network_connections[subnetwork_id]
-                allocation_model.subnetwork_demand[link] = zeros(n_demands)
-            end
+        add_secondary_network_demand!(allocation_model, p_independent)
+    else # Initialize subnetwork demands
+        n_priorities = length(p_independent.allocation.demand_priorities_all)
+
+        for link in p_independent.allocation.primary_network_connections[subnetwork_id]
+            allocation_model.secondary_network_demand[link] = zeros(n_priorities)
         end
     end