Add flow cost (#886)

Adding flow cost terms to linear objective functions. See also the added
documentation.

core/src/allocation.jl

@@ -661,6 +661,8 @@ function set_objective_priority!(
         ex = sum(problem[:F_abs])
     end
 
+    demand_max = 0.0
+
     for edge_id in edge_ids
         node_id_user = edge_id[2]
         if graph[node_id_user].type != :user
@@ -669,6 +671,7 @@ function set_objective_priority!(
 
         user_idx = findsorted(node_id, node_id_user)
         d = demand[user_idx][priority_idx](t)
+        demand_max = max(demand_max, d)
         F_edge = F[edge_id]
 
         if objective_type == :quadratic_absolute
@@ -707,6 +710,22 @@ function set_objective_priority!(
             error("Invalid allocation objective type $objective_type.")
         end
     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

core/src/config.jl

@@ -123,7 +123,7 @@ end
 @option struct Allocation <: TableOption
     timestep::Union{Float64, Nothing} = nothing
     use_allocation::Bool = false
-    objective_type::String = "quadratic_relative"
+    objective_type::String = "linear_absolute"
 end
 
 @option @addnodetypes struct Toml <: TableOption

core/test/allocation_test.jl

@@ -18,8 +18,16 @@
     allocation_model = p.allocation_models[1]
     Ribasim.allocate!(p, allocation_model, 0.0)
 
+    F = allocation_model.problem[:F]
+    @test JuMP.value(F[(NodeID(2), NodeID(6))]) ≈ 0.0
+    @test JuMP.value(F[(NodeID(2), NodeID(10))]) ≈ 0.5
+    @test JuMP.value(F[(NodeID(8), NodeID(12))]) ≈ 0.0
+    @test JuMP.value(F[(NodeID(6), NodeID(8))]) ≈ 0.0
+    @test JuMP.value(F[(NodeID(1), NodeID(2))]) ≈ 0.5
+    @test JuMP.value(F[(NodeID(6), NodeID(11))]) ≈ 0.0
+
     allocated = p.user.allocated
-    @test allocated[1] ≈ [0.0, 0.0]
+    @test allocated[1] ≈ [0.0, 0.5]
     @test allocated[2] ≈ [4.0, 0.0]
     @test allocated[3] ≈ [0.0, 0.0]
 end
@@ -74,8 +82,15 @@ end
     objective = JuMP.objective_function(problem)
     @test objective isa JuMP.AffExpr # Affine expression
     @test :F_abs in keys(problem.obj_dict)
+    F = problem[:F]
     F_abs = problem[:F_abs]
-    @test objective == F_abs[NodeID(5)] + F_abs[NodeID(6)]
+
+    @test objective.terms[F_abs[NodeID(5)]] == 1.0
+    @test objective.terms[F_abs[NodeID(6)]] == 1.0
+    @test objective.terms[F[(NodeID(4), NodeID(6))]] ≈ 0.125
+    @test objective.terms[F[(NodeID(1), NodeID(2))]] ≈ 0.125
+    @test objective.terms[F[(NodeID(4), NodeID(5))]] ≈ 0.125
+    @test objective.terms[F[(NodeID(2), NodeID(4))]] ≈ 0.125
 
     config = Ribasim.Config(toml_path; allocation_objective_type = "linear_relative")
     model = Ribasim.run(config)
@@ -84,6 +99,13 @@ end
     objective = JuMP.objective_function(problem)
     @test objective isa JuMP.AffExpr # Affine expression
     @test :F_abs in keys(problem.obj_dict)
+    F = problem[:F]
     F_abs = problem[:F_abs]
-    @test objective == F_abs[NodeID(5)] + F_abs[NodeID(6)]
+
+    @test objective.terms[F_abs[NodeID(5)]] == 1.0
+    @test objective.terms[F_abs[NodeID(6)]] == 1.0
+    @test objective.terms[F[(NodeID(4), NodeID(6))]] ≈ 62.585499316005475
+    @test objective.terms[F[(NodeID(1), NodeID(2))]] ≈ 62.585499316005475
+    @test objective.terms[F[(NodeID(4), NodeID(5))]] ≈ 62.585499316005475
+    @test objective.terms[F[(NodeID(2), NodeID(4))]] ≈ 62.585499316005475
 end

core/test/docs.toml

@@ -14,9 +14,9 @@ results_dir = "results" # required
 time = "basin/time.arrow"
 
 [allocation]
-timestep = 86400
-use_allocation = true
-objective_type = "quadratic_relative"
+timestep = 86400                   # optional (required if use_allocation = true), default 86400
+use_allocation = false             # optional, default false
+objective_type = "linear_absolute" # optional, default "linear_absolute"
 
 [solver]
 algorithm = "QNDF"  # optional, default "QNDF"

docs/core/allocation.qmd

@@ -112,19 +112,21 @@ $$
 $$
     \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left( 1 - \frac{F_{ij}}{d_j^p(t)}\right)^2
 $$
-- `linear_absolute`:
+- `linear_absolute` (default):
 $$
-    \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left| F_{ij} - d_j^p(t)\right|
+    \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
 $$
 - `linear_relative`:
 $$
-    \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left|1 - \frac{F_{ij}}{d_j^p(t)}\right|
+    \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
 $$
 
 To avoid division by $0$ errors, if a `*_relative` objective is used and a demand is $0$, the coefficient of the flow $F_{ij}$ is set to $0$.
 
 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.
 :::

docs/schema/Allocation.schema.json

@@ -25,7 +25,7 @@
         "objective_type": {
             "format": "default",
             "type": "string",
-            "default": "quadratic_relative"
+            "default": "linear_absolute"
         }
     },
     "required": [

docs/schema/Toml.schema.json

@@ -31,7 +31,7 @@
             "default": {
                 "timestep": null,
                 "use_allocation": false,
-                "objective_type": "quadratic_relative"
+                "objective_type": "linear_absolute"
             }
         },
         "solver": {

python/ribasim_testmodels/ribasim_testmodels/allocation.py

@@ -217,8 +217,8 @@ def subnetwork_model():
         time=pd.DataFrame(
             data={
                 "node_id": 1,
-                "flow_rate": np.arange(10, 0, -1),
-                "time": pd.to_datetime([f"2020-{i}-1 00:00:00" for i in range(1, 11)]),
+                "flow_rate": np.arange(10, 0, -2),
+                "time": pd.to_datetime([f"2020-{i}-1 00:00:00" for i in range(1, 6)]),
             },
         )
     )