From 45ad647dafad9a56fe8b29d90c123b180c52150e Mon Sep 17 00:00:00 2001
From: Jingru923 <47444880+Jingru923@users.noreply.github.com>
Date: Tue, 13 Feb 2024 14:27:43 +0100
Subject: [PATCH] change leq to eq constrain in flow conservation (#1083)

Fixes #954

changed the less than or equal constraints to equal constraints. I also
delete the flow cost term in the `linear_absolute` and `linear_relative`
objective functions.

This result in a change in the optimization solutions.
---
 core/src/allocation.jl       | 27 +++----------------
 core/test/allocation_test.jl | 51 ++++++++++++++++++++++++++----------
 docs/core/allocation.qmd     |  9 +++----
 3 files changed, 44 insertions(+), 43 deletions(-)

diff --git a/core/src/allocation.jl b/core/src/allocation.jl
index 5925b6677..82018d55e 100644
--- a/core/src/allocation.jl
+++ b/core/src/allocation.jl
@@ -495,14 +495,14 @@ Add the flow conservation constraints to the allocation problem.
 The constraint indices are user node IDs.
 
 Constraint:
-sum(flows out of node node) <= flows into node + flow from storage and vertical fluxes
+sum(flows out of node node) == flows into node + flow from storage and vertical fluxes
 """
 function add_constraints_flow_conservation!(
     problem::JuMP.Model,
     p::Parameters,
     allocation_network_id::Int,
 )::Nothing
-    (; graph, allocation) = p
+    (; graph) = p
     F = problem[:F]
     node_ids = graph[].node_ids[allocation_network_id]
     node_ids_conservation =
@@ -518,7 +518,7 @@ function add_constraints_flow_conservation!(
         sum([
             F[(node_id, outneighbor_id)] for
             outneighbor_id in outflow_ids_allocation(graph, node_id)
-        ]) <= sum([
+        ]) == sum([
             F[(inneighbor_id, node_id)] for
             inneighbor_id in inflow_ids_allocation(graph, node_id)
         ]),
@@ -824,14 +824,11 @@ function set_objective_priority!(
         ex = sum(problem[:F_abs])
     end
 
-    demand_max = 0.0
-
     # Terms for subnetworks as users
     if is_main_network(allocation_network_id)
         for connections_subnetwork in main_network_connections
             for connection in connections_subnetwork
                 d = subnetwork_demands[connection][priority_idx]
-                demand_max = max(demand_max, d)
                 add_user_term!(ex, connection, objective_type, d, allocation_model)
             end
         end
@@ -853,25 +850,9 @@ function set_objective_priority!(
             d = get_user_demand(user, node_id_user, priority_idx)
         end
 
-        demand_max = max(demand_max, d)
         add_user_term!(ex, edge_id, objective_type, d, allocation_model)
     end
 
-    # Add flow cost
-    if objective_type == :linear_absolute
-        cost_per_flow = 0.5 / length(F)
-        for flow in F
-            JuMP.add_to_expression!(ex, cost_per_flow * flow)
-        end
-    elseif objective_type == :linear_relative
-        if demand_max > 0.0
-            cost_per_flow = 0.5 / (demand_max * length(F))
-            for flow in F
-                JuMP.add_to_expression!(ex, cost_per_flow * flow)
-            end
-        end
-    end
-
     new_objective = JuMP.@expression(problem, ex)
     JuMP.@objective(problem, Min, new_objective)
     return nothing
@@ -1128,7 +1109,7 @@ function allocate!(
         @debug JuMP.solution_summary(problem)
         if JuMP.termination_status(problem) !== JuMP.OPTIMAL
             (; allocation_network_id) = allocation_model
-            priority = priorities[priority_index]
+            priority = priorities[priority_idx]
             error(
                 "Allocation of subnetwork $allocation_network_id, priority $priority coudn't find optimal solution.",
             )
diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl
index 63d504143..c39eb30e8 100644
--- a/core/test/allocation_test.jl
+++ b/core/test/allocation_test.jl
@@ -47,7 +47,7 @@
     @test Ribasim.get_user_demand(user, NodeID(:User, 11), 2) ≈ π
 end
 
-@testitem "Allocation objective types" begin
+@testitem "Allocation objective: quadratic absolute" skip = true begin
     using DataFrames: DataFrame
     using SciMLBase: successful_retcode
     using Ribasim: NodeID
@@ -72,6 +72,17 @@ end
         F[(NodeID(:Basin, 4), NodeID(:User, 6))],
         F[(NodeID(:Basin, 4), NodeID(:User, 6))],
     ) in keys(objective.terms) # F[4,6]^2 term
+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")
     model = Ribasim.run(config)
@@ -89,6 +100,17 @@ end
         F[(NodeID(:Basin, 4), NodeID(:User, 6))],
         F[(NodeID(:Basin, 4), NodeID(:User, 6))],
     ) in keys(objective.terms) # F[4,6]^2 term
+end
+
+@testitem "Allocation objective: linear 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 = "linear_absolute")
     model = Ribasim.run(config)
@@ -102,10 +124,17 @@ end
 
     @test objective.terms[F_abs[NodeID(:User, 5)]] == 1.0
     @test objective.terms[F_abs[NodeID(:User, 6)]] == 1.0
-    @test objective.terms[F[(NodeID(:Basin, 4), NodeID(:User, 6))]] ≈ 0.125
-    @test objective.terms[F[(NodeID(:FlowBoundary, 1), NodeID(:Basin, 2))]] ≈ 0.125
-    @test objective.terms[F[(NodeID(:Basin, 4), NodeID(:User, 5))]] ≈ 0.125
-    @test objective.terms[F[(NodeID(:Basin, 2), NodeID(:Basin, 4))]] ≈ 0.125
+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")
     model = Ribasim.run(config)
@@ -119,11 +148,6 @@ end
 
     @test objective.terms[F_abs[NodeID(:User, 5)]] == 1.0
     @test objective.terms[F_abs[NodeID(:User, 6)]] == 1.0
-    @test objective.terms[F[(NodeID(:Basin, 4), NodeID(:User, 6))]] ≈ 62.585499316005475
-    @test objective.terms[F[(NodeID(:FlowBoundary, 1), NodeID(:Basin, 2))]] ≈
-          62.585499316005475
-    @test objective.terms[F[(NodeID(:Basin, 4), NodeID(:User, 5))]] ≈ 62.585499316005475
-    @test objective.terms[F[(NodeID(:Basin, 2), NodeID(:Basin, 4))]] ≈ 62.585499316005475
 end
 
 @testitem "Allocation with controlled fractional flow" begin
@@ -259,8 +283,7 @@ end
     end
 
     @test subnetwork_demands[(NodeID(:Basin, 2), NodeID(:Pump, 11))] ≈ [4.0, 4.0, 0.0]
-    @test subnetwork_demands[(NodeID(:Basin, 6), NodeID(:Pump, 24))] ≈
-          [0.001333333333, 0.0, 0.0]
+    @test subnetwork_demands[(NodeID(:Basin, 6), NodeID(:Pump, 24))] ≈ [0.004, 0.0, 0.0]
     @test subnetwork_demands[(NodeID(:Basin, 10), NodeID(:Pump, 38))] ≈
           [0.001, 0.002, 0.002]
 
@@ -283,9 +306,9 @@ end
     Ribasim.update_allocation!((; p, t))
 
     @test subnetwork_allocateds[NodeID(:Basin, 2), NodeID(:Pump, 11)] ≈
-          [4.0, 0.49766666, 0.0]
+          [4.0, 0.49500000, 0.0]
     @test subnetwork_allocateds[NodeID(:Basin, 6), NodeID(:Pump, 24)] ≈
-          [0.00133333333, 0.0, 0.0]
+          [0.00399999999, 0.0, 0.0]
     @test subnetwork_allocateds[NodeID(:Basin, 10), NodeID(:Pump, 38)] ≈ [0.001, 0.0, 0.0]
 
     @test user.allocated[2] ≈ [4.0, 0.0, 0.0]
diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd
index 6bf3d266d..8c235b5e0 100644
--- a/docs/core/allocation.qmd
+++ b/docs/core/allocation.qmd
@@ -130,11 +130,11 @@ $$
 $$
 - `linear_absolute` (default):
 $$
-    \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left| F_{ij} - d_j^p(t)\right| + c \sum_{e \in E_S} F_e
+    \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left| F_{ij} - d_j^p(t)\right|
 $$
 - `linear_relative`:
 $$
-    \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left|1 - \frac{F_{ij}}{d_j^p(t)}\right| + c \sum_{e \in E_S} F_e
+    \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left|1 - \frac{F_{ij}}{d_j^p(t)}\right|
 $$
 
 :::{.callout-note}
@@ -146,8 +146,6 @@ To avoid division by $0$ errors, if a `*_relative` objective is used and a deman
 
 For `*_absolute` objectives the optimizer cares about the actual amount of water allocated to a user, for `*_relative` objectives it cares about the fraction of the demand allocated to the user. For `quadratic_*` objectives the optimizer cares about avoiding large shortages, for `linear_*` objectives it treats all deviations equally.
 
-The second sum in the `linear_*` objectives adds a small cost to using flows. This incentivizes the solver to use as little flow as possible. The cost $c > 0$ is small enough such that it is always better to bring water to users than to not use flow at all. This can be achieved for `linear_*` objectives but not for `quadratic_*` objectives, and therefore this cost term is only added to the former. Therefore the `linear_*` objectives make the solver more conservative with flow than the `quadratic_*` objectives.
-
 :::{.callout-note}
 These options for objectives for allocation to users have not been tested thoroughly, and might change in the future.
 :::
@@ -161,9 +159,8 @@ In the future new optimization objectives will be introduced, for demands of bas
 ## The optimization constraints
 - Flow conservation: For the basins in the allocation graph we have that
 $$
-    \sum_{j=1}^{n'} F_{kj} \le \sum_{i=1}^{n'} F_{ik}, \quad \forall k \in B_S.
+    \sum_{j=1}^{n'} F_{kj} = \sum_{i=1}^{n'} F_{ik}, \quad \forall k \in B_S.
 $$  {#eq-flowconservationconstraint}
-Note that we do not require equality here; in the allocation we do not mind that excess flow is 'forgotten' if it cannot contribute to the allocation to the users.
 - Capacity: the flows over the edges are positive and bounded by the edge capacity:
 $$
     F_{ij} \le \left(C_S\right)_{ij}, \quad \forall(i,j) \in E_S.