Make LevelDemand min_level and max_level optional (#1430)

As proposed in
#1428 (comment)

This makes it easier to use LevelDemand in case you only care about a
shortage or surplus situation.

---------

Co-authored-by: Bart de Koning <[email protected]>

core/src/read.jl

@@ -93,7 +93,7 @@ function parse_static_and_time(
 
     errors = false
     t_end = seconds_since(config.endtime, config.starttime)
-    trivial_timespan = [nextfloat(-Inf), prevfloat(Inf)]
+    trivial_timespan = [0.0, prevfloat(Inf)]
 
     for (node_idx, node_id) in enumerate(node_ids)
         if node_id in static_node_ids
@@ -820,7 +820,7 @@ function UserDemand(db::DB, config::Config)::UserDemand
     realized_bmi = zeros(n_user)
     demand = zeros(n_user, n_priority)
     demand_reduced = zeros(n_user, n_priority)
-    trivial_timespan = [nextfloat(-Inf), prevfloat(Inf)]
+    trivial_timespan = [0.0, prevfloat(Inf)]
     demand_itp = [
         [LinearInterpolation(zeros(2), trivial_timespan) for i in eachindex(priorities)] for j in eachindex(node_ids)
     ]
@@ -885,6 +885,7 @@ function LevelDemand(db::DB, config::Config)::LevelDemand
         static,
         time,
         time_interpolatables = [:min_level, :max_level],
+        defaults = (; min_level = -Inf, max_level = Inf),
     )
 
     if !valid

core/src/schema.jl

@@ -250,16 +250,16 @@ end
 
 @version LevelDemandStaticV1 begin
     node_id::Int32
-    min_level::Float64
-    max_level::Float64
+    min_level::Union{Missing, Float64}
+    max_level::Union{Missing, Float64}
     priority::Int32
 end
 
 @version LevelDemandTimeV1 begin
     node_id::Int32
     time::DateTime
-    min_level::Float64
-    max_level::Float64
+    min_level::Union{Missing, Float64}
+    max_level::Union{Missing, Float64}
     priority::Int32
 end
 

core/src/util.jl

@@ -31,17 +31,10 @@ function get_storage_from_level(basin::Basin, state_idx::Int, level::Float64)::F
     storage_discrete = basin.storage[state_idx]
     area_discrete = basin.area[state_idx]
     level_discrete = basin.level[state_idx]
-    bottom = first(level_discrete)
-
-    if level < bottom
-        node_id = basin.node_id.values[state_idx]
-        @error "The level $level of $node_id is lower than the bottom of this basin; $bottom."
-        return NaN
-    end
 
     level_lower_index = searchsortedlast(level_discrete, level)
 
-    # If the level is equal to the bottom then the storage is 0
+    # If the level is at or below the bottom then the storage is 0
     if level_lower_index == 0
         return 0.0
     end
@@ -77,7 +70,10 @@ function get_storages_from_levels(basin::Basin, levels::Vector)::Vector{Float64}
 
     for (i, level) in enumerate(levels)
         storage = get_storage_from_level(basin, i, level)
-        if isnan(storage)
+        bottom = first(basin.level[i])
+        node_id = basin.node_id.values[i]
+        if level < bottom
+            @error "The initial level ($level) of $node_id is below the bottom ($bottom)."
             errors = true
         end
         storages[i] = storage

core/test/allocation_test.jl

@@ -313,6 +313,8 @@ end
 end
 
 @testitem "Allocation level control" begin
+    import JuMP
+
     toml_path = normpath(@__DIR__, "../../generated_testmodels/level_demand/ribasim.toml")
     @test ispath(toml_path)
     model = Ribasim.run(toml_path)
@@ -377,6 +379,14 @@ end
     stage_6_start_idx = findfirst(stage_6)
     u_stage_6(τ) = storage[stage_6_start_idx]
     @test storage[stage_6] ≈ u_stage_6.(t[stage_6]) rtol = 1e-4
+
+    # Isolated LevelDemand + Basin pair to test optional min_level
+    problem = allocation.allocation_models[2].problem
+    @test JuMP.value(only(problem[:F_basin_in])) == 0.0
+    @test JuMP.value(only(problem[:F_basin_out])) == 0.0
+    q = JuMP.normalized_rhs(only(problem[:basin_outflow]))
+    storage_surplus = 1000.0  # Basin #7 is 1000 m2 and 1 m above LevelDemand max_level
+    @test q ≈ storage_surplus / Δt_allocation
 end
 
 @testitem "Flow demand" begin

core/test/utils_test.jl

@@ -114,14 +114,17 @@ end
     @test length(logger.logs) == 1
     @test logger.logs[1].level == Error
     @test logger.logs[1].message ==
-          "The level -1.0 of Basin #1 is lower than the bottom of this basin; 0.0."
+          "The initial level (-1.0) of Basin #1 is below the bottom (0.0)."
 
     # Converting from storages to levels and back should return the same storages
     storages = range(0.0, 2 * storage[end], 50)
     levels = [Ribasim.get_area_and_level(basin, 1, s)[2] for s in storages]
     storages_ = [Ribasim.get_storage_from_level(basin, 1, l) for l in levels]
-
     @test storages ≈ storages_
+
+    # At or below bottom the storage is 0
+    @test Ribasim.get_storage_from_level(basin, 1, 0.0) == 0.0
+    @test Ribasim.get_storage_from_level(basin, 1, -1.0) == 0.0
 end
 
 @testitem "Expand logic_mapping" begin

docs/core/usage.qmd

@@ -479,14 +479,17 @@ min_level     | Float64  | $m$          | -
 
 A `LevelDemand` node associates a minimum and a maximum level with connected basins to be used by the allocation algorithm.
 Below the minimum level the basin has a demand of the supplied priority,
-above the maximum level the basin acts as a source, used by all nodes with demands in order of priority.
-The same `LevelDemand` node can be used for basins in different subnetworks.
+above the maximum level the basin has a surplus and acts as a source, used by all nodes with demands in order of priority.
+The same `LevelDemand` node can be used for Basins in different subnetworks.
+
+Both `min_level` and `max_level` are optional, to be able to handle only the demand or surplus side.
+If both are missing, `LevelDemand` won't have any effects on allocation.
 
 column        | type    | unit         | restriction
 ------------- | ------- | ------------ | -----------
 node_id       | Int32   | -            | sorted
-min_level     | Float64 | $m$          | -
-max_level     | Float64 | $m$          | -
+min_level     | Float64 | $m$          | (optional, default -Inf)
+max_level     | Float64 | $m$          | (optional, default Inf)
 priority      | Int32   | -            | positive
 
 ## LevelDemand / time

python/ribasim/ribasim/schemas.py

@@ -115,16 +115,16 @@ class LevelBoundaryTimeSchema(_BaseSchema):
 
 class LevelDemandStaticSchema(_BaseSchema):
     node_id: Series[Int32] = pa.Field(nullable=False, default=0)
-    min_level: Series[float] = pa.Field(nullable=False)
-    max_level: Series[float] = pa.Field(nullable=False)
+    min_level: Series[float] = pa.Field(nullable=True)
+    max_level: Series[float] = pa.Field(nullable=True)
     priority: Series[Int32] = pa.Field(nullable=False, default=0)
 
 
 class LevelDemandTimeSchema(_BaseSchema):
     node_id: Series[Int32] = pa.Field(nullable=False, default=0)
     time: Series[Timestamp] = pa.Field(nullable=False)
-    min_level: Series[float] = pa.Field(nullable=False)
-    max_level: Series[float] = pa.Field(nullable=False)
+    min_level: Series[float] = pa.Field(nullable=True)
+    max_level: Series[float] = pa.Field(nullable=True)
     priority: Series[Int32] = pa.Field(nullable=False, default=0)
 
 

python/ribasim_testmodels/ribasim_testmodels/allocation.py

@@ -883,7 +883,7 @@ def subnetworks_with_sources_model() -> Model:
 
 
 def level_demand_model() -> Model:
-    """Small model with a LevelDemand."""
+    """Small model with LevelDemand nodes."""
 
     model = Model(
         starttime="2020-01-01",
@@ -922,12 +922,24 @@ def level_demand_model() -> Model:
         [basin.Profile(area=1000.0, level=[0.0, 1.0]), basin.State(level=[0.5])],
     )
 
+    # Isolated LevelDemand + Basin pair to test optional min_level
+    model.level_demand.add(
+        Node(6, Point(3, -1), subnetwork_id=3),
+        [level_demand.Static(max_level=[1.0], priority=1)],
+    )
+    model.basin.add(
+        Node(7, Point(3, 0), subnetwork_id=3),
+        [basin.Profile(area=1000.0, level=[0.0, 1.0]), basin.State(level=[2.0])],
+    )
+
     model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2)
     model.edge.add(model.basin[2], model.user_demand[3])
     model.edge.add(model.level_demand[4], model.basin[2])
     model.edge.add(model.user_demand[3], model.basin[5])
     model.edge.add(model.level_demand[4], model.basin[5])
 
+    model.edge.add(model.level_demand[6], model.basin[7])
+
     return model