From 38c981c144689a699ae576e9f536f1299a548e20 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 6 Feb 2024 11:50:06 +0100 Subject: [PATCH 01/37] Add node AllocationLevelControl in Julia core --- core/src/create.jl | 29 +++++++++++++++++++++++++++++ core/src/solve.jl | 8 ++++++++ core/src/validation.jl | 15 +++++++++++++++ 3 files changed, 52 insertions(+) diff --git a/core/src/create.jl b/core/src/create.jl index 0c571e7af..4858856b7 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -758,6 +758,33 @@ function User(db::DB, config::Config)::User ) end +function AllocationLevelControl(db::DB, config::Config)::AllocationLevelControl + static = load_structvector(db, config, AllocationLevelControlStaticV1) + time = load_structvector(db, config, AllocationLevelControlTimeV1) + + parsed_parameters, valid = parse_static_and_time( + db, + config, + "AllocationLevelControl"; + static, + time, + time_interpolatables = [:target_level], + ) + + if !valid + error("Errors occurred when parsing AllocationLevelControl data.") + end + + target_level = zeros(length(parsed_parameters.node_id)) + + return AllocationLevelControl( + NodeID.(parsed_parameters.node_id), + target_level, + parsed_parameters.target_level, + parsed_parameters.priority, + ) +end + function Subgrid(db::DB, config::Config, basin::Basin)::Subgrid node_to_basin = Dict(node_id => index for (index, node_id) in enumerate(basin.node_id)) tables = load_structvector(db, config, BasinSubgridV1) @@ -842,6 +869,7 @@ function Parameters(db::DB, config::Config)::Parameters discrete_control = DiscreteControl(db, config) pid_control = PidControl(db, config, chunk_sizes) user = User(db, config) + allocation_level_control = AllocationLevelControl(db, config) basin = Basin(db, config, chunk_sizes) subgrid_level = Subgrid(db, config, basin) @@ -875,6 +903,7 @@ function Parameters(db::DB, config::Config)::Parameters discrete_control, pid_control, user, + allocation_level_control, Dict{Int, Symbol}(), subgrid_level, ) diff --git a/core/src/solve.jl b/core/src/solve.jl index 317453640..eec444f4a 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -474,6 +474,13 @@ struct User <: AbstractParameterNode end end +struct AllocationLevelControl + node_id::Vector{NodeID} + target_level::Vector{Float64} + target_level_itp::Vector{LinearInterpolation} + priority::Vector{Int} +end + "Subgrid linearly interpolates basin levels." struct Subgrid basin_index::Vector{Int} @@ -516,6 +523,7 @@ struct Parameters{T, C1, C2} discrete_control::DiscreteControl pid_control::PidControl{T} user::User + allocation_level_control::AllocationLevelControl lookup::Dict{Int, Symbol} subgrid::Subgrid end diff --git a/core/src/validation.jl b/core/src/validation.jl index adfcc1da6..fb61fca40 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -25,6 +25,8 @@ @schema "ribasim.outlet.static" OutletStatic @schema "ribasim.user.static" UserStatic @schema "ribasim.user.time" UserTime +@schema "ribasim.allocationlevelcontrol.static" AllocationLevelControlStatic +@schema "ribasim.allocationlevelcontrol.time" AllocationLevelControlTime const delimiter = " / " tablename(sv::Type{SchemaVersion{T, N}}) where {T, N} = tablename(sv()) @@ -330,6 +332,19 @@ end priority::Int end +@version AllocationLevelControlStaticV1 begin + node_id::Int + target_level::Float64 + priority::Int +end + +@version AllocationLevelControlTimeV1 begin + node_id::Int + time::DateTime + target_level::Float64 + priority::Int +end + function variable_names(s::Any) filter(x -> !(x in (:node_id, :control_state)), fieldnames(s)) end From fdc6ccb10853500a7c0fc31f838540fdcaaf9762 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 6 Feb 2024 13:10:36 +0100 Subject: [PATCH 02/37] Add AllocationLevelControl node to QGIS plugin --- ribasim_qgis/core/nodes.py | 42 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/ribasim_qgis/core/nodes.py b/ribasim_qgis/core/nodes.py index 84f240139..1c034e857 100644 --- a/ribasim_qgis/core/nodes.py +++ b/ribasim_qgis/core/nodes.py @@ -215,6 +215,11 @@ def renderer(self) -> QgsCategorizedSymbolRenderer: "DiscreteControl": (QColor("black"), "DiscreteControl", shape.Star), "PidControl": (QColor("black"), "PidControl", shape.Cross2), "User": (QColor("green"), "User", shape.Square), + "AllocationLevelControl": ( + QColor("black"), + "AllocationLevelControl", + shape.Circle, + ), # All other nodes, or incomplete input "": (QColor("white"), "", shape.Circle), } @@ -780,6 +785,43 @@ def attributes(cls) -> list[QgsField]: ] +class AllocationLevelControl(Input): + @classmethod + def input_type(cls) -> str: + return "AllocationLevelControl / static" + + @classmethod + def geometry_type(cls) -> str: + return "No Geometry" + + @classmethod + def attributes(cls) -> list[QgsField]: + return [ + QgsField("node_id", QVariant.Int), + QgsField("target_level", QVariant.Double), + QgsField("priority", QVariant.Int), + ] + + +class AllocationLevelControlTime(Input): + @classmethod + def input_type(cls) -> str: + return "AllocationLevelControl / time" + + @classmethod + def geometry_type(cls) -> str: + return "No Geometry" + + @classmethod + def attributes(cls) -> list[QgsField]: + return [ + QgsField("node_id", QVariant.Int), + QgsField("time", QVariant.DateTime), + QgsField("target_level", QVariant.Double), + QgsField("priority", QVariant.Int), + ] + + NODES: dict[str, type[Input]] = { cls.input_type(): cls # type: ignore[type-abstract] # mypy doesn't see that all classes are concrete. for cls in Input.__subclasses__() From 31d660321fd020dcbaed06dcc8c93d6298ca7942 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 6 Feb 2024 13:17:11 +0100 Subject: [PATCH 03/37] Add validation rules for AllocationLevelControl --- core/src/validation.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/core/src/validation.jl b/core/src/validation.jl index fb61fca40..d04dabb8f 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -65,6 +65,7 @@ neighbortypes(nodetype::Symbol) = neighbortypes(Val(nodetype)) neighbortypes(::Val{:pump}) = Set((:basin, :fractional_flow, :terminal, :level_boundary)) neighbortypes(::Val{:outlet}) = Set((:basin, :fractional_flow, :terminal, :level_boundary)) neighbortypes(::Val{:user}) = Set((:basin, :fractional_flow, :terminal, :level_boundary)) +neighbortypes(::Val{:AllocationLevelControl}) = Set((:basin,)) neighbortypes(::Val{:basin}) = Set(( :linear_resistance, :tabulated_rating_curve, @@ -119,6 +120,7 @@ n_neighbor_bounds_flow(::Val{:Terminal}) = n_neighbor_bounds(1, typemax(Int), 0, n_neighbor_bounds_flow(::Val{:PidControl}) = n_neighbor_bounds(0, 0, 0, 0) n_neighbor_bounds_flow(::Val{:DiscreteControl}) = n_neighbor_bounds(0, 0, 0, 0) n_neighbor_bounds_flow(::Val{:User}) = n_neighbor_bounds(1, 1, 1, 1) +n_neighbor_bounds_flow(::Val{:AllocationLevelControl}) = n_neighbor_bounds(0, 0, 0, 0) n_neighbor_bounds_flow(nodetype) = error("'n_neighbor_bounds_flow' not defined for $nodetype.") @@ -137,6 +139,7 @@ n_neighbor_bounds_control(::Val{:PidControl}) = n_neighbor_bounds(0, 1, 1, 1) n_neighbor_bounds_control(::Val{:DiscreteControl}) = n_neighbor_bounds(0, 0, 1, typemax(Int)) n_neighbor_bounds_control(::Val{:User}) = n_neighbor_bounds(0, 0, 0, 0) +n_neighbor_bounds_control(::Val{:AllocationLevelControl}) = n_neighbor_bounds(0, 0, 1, 1) n_neighbor_bounds_control(nodetype) = error("'n_neighbor_bounds_control' not defined for $nodetype.") From 06db1ac4e3b14b3da665fad40d570aed3d8d5eea Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 6 Feb 2024 13:35:31 +0100 Subject: [PATCH 04/37] add AllocationLevelControl to Ribasim python --- .../AllocationLevelControlStatic.schema.json | 32 ++++++++++++++++ .../AllocationLevelControlTime.schema.json | 37 +++++++++++++++++++ docs/schema/Toml.schema.json | 8 ++++ .../allocation_level_control.schema.json | 35 ++++++++++++++++++ docs/schema/root.schema.json | 6 +++ python/ribasim/ribasim/__init__.py | 2 + python/ribasim/ribasim/config.py | 13 +++++++ python/ribasim/ribasim/models.py | 17 +++++++++ .../ribasim_testmodels/allocation.py | 4 ++ 9 files changed, 154 insertions(+) create mode 100644 docs/schema/AllocationLevelControlStatic.schema.json create mode 100644 docs/schema/AllocationLevelControlTime.schema.json create mode 100644 docs/schema/allocation_level_control.schema.json diff --git a/docs/schema/AllocationLevelControlStatic.schema.json b/docs/schema/AllocationLevelControlStatic.schema.json new file mode 100644 index 000000000..b04fc69a7 --- /dev/null +++ b/docs/schema/AllocationLevelControlStatic.schema.json @@ -0,0 +1,32 @@ +{ + "$schema": "https://json-schema.org/draft/2020-12/schema", + "$id": "https://deltares.github.io/Ribasim/schema/AllocationLevelControlStatic.schema.json", + "title": "AllocationLevelControlStatic", + "description": "A AllocationLevelControlStatic object based on Ribasim.AllocationLevelControlStaticV1", + "type": "object", + "properties": { + "node_id": { + "format": "default", + "type": "integer" + }, + "target_level": { + "format": "double", + "type": "number" + }, + "priority": { + "format": "default", + "type": "integer" + }, + "remarks": { + "description": "a hack for pandera", + "type": "string", + "format": "default", + "default": "" + } + }, + "required": [ + "node_id", + "target_level", + "priority" + ] +} diff --git a/docs/schema/AllocationLevelControlTime.schema.json b/docs/schema/AllocationLevelControlTime.schema.json new file mode 100644 index 000000000..5478922de --- /dev/null +++ b/docs/schema/AllocationLevelControlTime.schema.json @@ -0,0 +1,37 @@ +{ + "$schema": "https://json-schema.org/draft/2020-12/schema", + "$id": "https://deltares.github.io/Ribasim/schema/AllocationLevelControlTime.schema.json", + "title": "AllocationLevelControlTime", + "description": "A AllocationLevelControlTime object based on Ribasim.AllocationLevelControlTimeV1", + "type": "object", + "properties": { + "node_id": { + "format": "default", + "type": "integer" + }, + "time": { + "format": "date-time", + "type": "string" + }, + "target_level": { + "format": "double", + "type": "number" + }, + "priority": { + "format": "default", + "type": "integer" + }, + "remarks": { + "description": "a hack for pandera", + "type": "string", + "format": "default", + "default": "" + } + }, + "required": [ + "node_id", + "time", + "target_level", + "priority" + ] +} diff --git a/docs/schema/Toml.schema.json b/docs/schema/Toml.schema.json index 1fce57eba..af232ee5a 100644 --- a/docs/schema/Toml.schema.json +++ b/docs/schema/Toml.schema.json @@ -76,6 +76,13 @@ "static": null } }, + "allocation_level_control": { + "$ref": "https://deltares.github.io/Ribasim/schema/allocation_level_control.schema.json", + "default": { + "static": null, + "time": null + } + }, "pid_control": { "$ref": "https://deltares.github.io/Ribasim/schema/pid_control.schema.json", "default": { @@ -170,6 +177,7 @@ "logging", "results", "terminal", + "allocation_level_control", "pid_control", "level_boundary", "pump", diff --git a/docs/schema/allocation_level_control.schema.json b/docs/schema/allocation_level_control.schema.json new file mode 100644 index 000000000..9edae20da --- /dev/null +++ b/docs/schema/allocation_level_control.schema.json @@ -0,0 +1,35 @@ +{ + "$schema": "https://json-schema.org/draft/2020-12/schema", + "$id": "https://deltares.github.io/Ribasim/schema/allocation_level_control.schema.json", + "title": "allocation_level_control", + "description": "A allocation_level_control object based on Ribasim.config.allocation_level_control", + "type": "object", + "properties": { + "static": { + "format": "default", + "anyOf": [ + { + "type": "null" + }, + { + "type": "string" + } + ], + "default": null + }, + "time": { + "format": "default", + "anyOf": [ + { + "type": "null" + }, + { + "type": "string" + } + ], + "default": null + } + }, + "required": [ + ] +} diff --git a/docs/schema/root.schema.json b/docs/schema/root.schema.json index e9089cc68..8bb74eb6a 100644 --- a/docs/schema/root.schema.json +++ b/docs/schema/root.schema.json @@ -1,6 +1,12 @@ { "$schema": "https://json-schema.org/draft/2020-12/schema", "properties": { + "allocationlevelcontrolstatic": { + "$ref": "allocationlevelcontrolstatic.schema.json" + }, + "allocationlevelcontroltime": { + "$ref": "allocationlevelcontroltime.schema.json" + }, "basinprofile": { "$ref": "basinprofile.schema.json" }, diff --git a/python/ribasim/ribasim/__init__.py b/python/ribasim/ribasim/__init__.py index e7b2a7de9..c853fb1ed 100644 --- a/python/ribasim/ribasim/__init__.py +++ b/python/ribasim/ribasim/__init__.py @@ -4,6 +4,7 @@ from ribasim import models, utils from ribasim.config import ( Allocation, + AllocationLevelControl, Basin, Compression, DiscreteControl, @@ -29,6 +30,7 @@ __all__ = [ "Allocation", + "AllocationLevelControl", "Basin", "DiscreteControl", "Compression", diff --git a/python/ribasim/ribasim/config.py b/python/ribasim/ribasim/config.py index 167f8fb64..01a69c88d 100644 --- a/python/ribasim/ribasim/config.py +++ b/python/ribasim/ribasim/config.py @@ -6,6 +6,8 @@ # These schemas are autogenerated from ribasim.schemas import ( # type: ignore + AllocationLevelControlStaticSchema, + AllocationLevelControlTimeSchema, BasinProfileSchema, BasinStateSchema, BasinStaticSchema, @@ -135,6 +137,17 @@ class User(NodeModel): ) +class AllocationLevelControl(NodeModel): + static: TableModel[AllocationLevelControlStaticSchema] = Field( + default_factory=TableModel[AllocationLevelControlStaticSchema], + json_schema_extra={"sort_keys": ["node_id", "priority"]}, + ) + time: TableModel[AllocationLevelControlTimeSchema] = Field( + default_factory=TableModel[AllocationLevelControlTimeSchema], + json_schema_extra={"sort_keys": ["node_id", "priority", "time"]}, + ) + + class FlowBoundary(NodeModel): static: TableModel[FlowBoundaryStaticSchema] = Field( default_factory=TableModel[FlowBoundaryStaticSchema], diff --git a/python/ribasim/ribasim/models.py b/python/ribasim/ribasim/models.py index d0e8c341b..3e7faa7bb 100644 --- a/python/ribasim/ribasim/models.py +++ b/python/ribasim/ribasim/models.py @@ -10,6 +10,21 @@ from ribasim.input_base import BaseModel +class AllocationLevelControlStatic(BaseModel): + node_id: int + target_level: float + priority: int + remarks: str = Field("", description="a hack for pandera") + + +class AllocationLevelControlTime(BaseModel): + node_id: int + time: datetime + target_level: float + priority: int + remarks: str = Field("", description="a hack for pandera") + + class BasinProfile(BaseModel): node_id: int area: float @@ -228,6 +243,8 @@ class UserTime(BaseModel): class Root(BaseModel): + allocationlevelcontrolstatic: AllocationLevelControlStatic | None = None + allocationlevelcontroltime: AllocationLevelControlTime | None = None basinprofile: BasinProfile | None = None basinstate: BasinState | None = None basinstatic: BasinStatic | None = None diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index a0bfa342d..eef04c6a5 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -1735,3 +1735,7 @@ def main_network_with_subnetworks_model(): ) return model + + +def allocation_level_control_model(): + return From 67f2cf203f3280ca0df35eb63e5443049133ea66 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 6 Feb 2024 15:18:39 +0100 Subject: [PATCH 05/37] add test model allocation_level_control --- core/src/allocation.jl | 3 +- core/src/validation.jl | 7 +- docs/core/usage.qmd | 2 +- python/ribasim/ribasim/geometry/node.py | 2 + python/ribasim/ribasim/model.py | 4 + .../ribasim_testmodels/__init__.py | 2 + .../ribasim_testmodels/allocation.py | 104 +++++++++++++++++- 7 files changed, 117 insertions(+), 7 deletions(-) diff --git a/core/src/allocation.jl b/core/src/allocation.jl index d49cdf0da..2c4b8e2dd 100644 --- a/core/src/allocation.jl +++ b/core/src/allocation.jl @@ -854,7 +854,8 @@ function set_objective_priority!( if objective_type in [:quadratic_absolute, :quadratic_relative] ex = JuMP.QuadExpr() elseif objective_type in [:linear_absolute, :linear_relative] - ex = sum(problem[:F_abs]) + ex = JuMP.AffExpr() + ex += sum(problem[:F_abs]) end demand_max = 0.0 diff --git a/core/src/validation.jl b/core/src/validation.jl index d04dabb8f..b3cbf26ab 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -65,7 +65,7 @@ neighbortypes(nodetype::Symbol) = neighbortypes(Val(nodetype)) neighbortypes(::Val{:pump}) = Set((:basin, :fractional_flow, :terminal, :level_boundary)) neighbortypes(::Val{:outlet}) = Set((:basin, :fractional_flow, :terminal, :level_boundary)) neighbortypes(::Val{:user}) = Set((:basin, :fractional_flow, :terminal, :level_boundary)) -neighbortypes(::Val{:AllocationLevelControl}) = Set((:basin,)) +neighbortypes(::Val{:allocation_level_control}) = Set((:basin,)) neighbortypes(::Val{:basin}) = Set(( :linear_resistance, :tabulated_rating_curve, @@ -125,7 +125,7 @@ n_neighbor_bounds_flow(nodetype) = error("'n_neighbor_bounds_flow' not defined for $nodetype.") n_neighbor_bounds_control(nodetype::Symbol) = n_neighbor_bounds_control(Val(nodetype)) -n_neighbor_bounds_control(::Val{:Basin}) = n_neighbor_bounds(0, 0, 0, typemax(Int)) +n_neighbor_bounds_control(::Val{:Basin}) = n_neighbor_bounds(0, 1, 0, typemax(Int)) n_neighbor_bounds_control(::Val{:LinearResistance}) = n_neighbor_bounds(0, 1, 0, 0) n_neighbor_bounds_control(::Val{:ManningResistance}) = n_neighbor_bounds(0, 1, 0, 0) n_neighbor_bounds_control(::Val{:TabulatedRatingCurve}) = n_neighbor_bounds(0, 1, 0, 0) @@ -139,7 +139,8 @@ n_neighbor_bounds_control(::Val{:PidControl}) = n_neighbor_bounds(0, 1, 1, 1) n_neighbor_bounds_control(::Val{:DiscreteControl}) = n_neighbor_bounds(0, 0, 1, typemax(Int)) n_neighbor_bounds_control(::Val{:User}) = n_neighbor_bounds(0, 0, 0, 0) -n_neighbor_bounds_control(::Val{:AllocationLevelControl}) = n_neighbor_bounds(0, 0, 1, 1) +n_neighbor_bounds_control(::Val{:AllocationLevelControl}) = + n_neighbor_bounds(0, 0, 1, typemax(Int)) n_neighbor_bounds_control(nodetype) = error("'n_neighbor_bounds_control' not defined for $nodetype.") diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index fe1cb5013..25c18fbf7 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -470,7 +470,7 @@ node_id | Int | - | sorted active | Bool | - | (optional, default true) demand | Float64 | $m^3 s^{-1}$ | - return_factor | Float64 | - | between [0 - 1] -min_level | Float64 | $m$ | (optional) +min_level | Float64 | $m$ | - priority | Int | - | sorted per node id ## User / time diff --git a/python/ribasim/ribasim/geometry/node.py b/python/ribasim/ribasim/geometry/node.py index 5c3752b00..2fc31ce66 100644 --- a/python/ribasim/ribasim/geometry/node.py +++ b/python/ribasim/ribasim/geometry/node.py @@ -207,6 +207,7 @@ def plot(self, ax=None, zorder=None) -> Any: "DiscreteControl": "*", "PidControl": "x", "User": "s", + "AllocationLevelControl": "o", "": "o", } @@ -224,6 +225,7 @@ def plot(self, ax=None, zorder=None) -> Any: "DiscreteControl": "k", "PidControl": "k", "User": "g", + "AllocationLevelControl": "k", "": "k", } assert self.df is not None diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py index 05057a58c..0d36610f5 100644 --- a/python/ribasim/ribasim/model.py +++ b/python/ribasim/ribasim/model.py @@ -19,6 +19,7 @@ from ribasim.config import ( Allocation, + AllocationLevelControl, Basin, DiscreteControl, FlowBoundary, @@ -171,6 +172,9 @@ class Model(FileModel): logging: Logging = Logging() allocation: Allocation = Field(default_factory=Allocation) + allocation_level_control: AllocationLevelControl = Field( + default_factory=AllocationLevelControl + ) basin: Basin = Field(default_factory=Basin) fractional_flow: FractionalFlow = Field(default_factory=FractionalFlow) level_boundary: LevelBoundary = Field(default_factory=LevelBoundary) diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index 4616d2746..42adc6c6f 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -7,6 +7,7 @@ import ribasim_testmodels from ribasim_testmodels.allocation import ( allocation_example_model, + allocation_level_control_model, fractional_flow_subnetwork_model, looped_subnetwork_model, main_network_with_subnetworks_model, @@ -54,6 +55,7 @@ __all__ = [ "allocation_example_model", + "allocation_level_control_model", "backwater_model", "basic_arrow_model", "basic_model", diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index eef04c6a5..38dc458b5 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -1081,7 +1081,9 @@ def allocation_example_model(): def main_network_with_subnetworks_model(): - # Set ip the nodes: + """Generate a model which consists of a main network and multiple connected subnetworks.""" + + # Set up the nodes: xy = np.array( [ (0.0, -1.0), @@ -1738,4 +1740,102 @@ def main_network_with_subnetworks_model(): def allocation_level_control_model(): - return + # Set up the nodes: + xy = np.array( + [ + (0.0, 0.0), # 1: FlowBoundary + (1.0, 0.0), # 2: Basin + (2.0, 0.0), # 3: User + (1.0, -1.0), # 4: AllocationLevelControl + ] + ) + node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1]) + + node_type = ["FlowBoundary", "Basin", "User", "AllocationLevelControl"] + + # Make sure the feature id starts at 1: explicitly give an index. + node = ribasim.Node( + df=gpd.GeoDataFrame( + data={ + "type": node_type, + "allocation_network_id": 4 * [1], + }, + index=pd.Index(np.arange(len(xy)) + 1, name="fid"), + geometry=node_xy, + crs="EPSG:28992", + ) + ) + + # Setup the edges: + from_id = np.array([1, 2, 4, 3]) + to_id = np.array([2, 3, 2, 2]) + edge_type = ["flow", "flow", "control", "flow"] + allocation_network_id = [1, None, None, None] + + lines = node.geometry_from_connectivity(from_id.tolist(), to_id.tolist()) + edge = ribasim.Edge( + df=gpd.GeoDataFrame( + data={ + "from_node_id": from_id, + "to_node_id": to_id, + "edge_type": edge_type, + "allocation_network_id": allocation_network_id, + }, + geometry=lines, + crs="EPSG:28992", + ) + ) + + # Setup basin + profile = pd.DataFrame(data={"node_id": 2, "area": 1e3, "level": [0.0, 1.0]}) + static = pd.DataFrame( + data={ + "node_id": [2], + "drainage": 0.0, + "potential_evaporation": 0.0, + "infiltration": 0.0, + "precipitation": 0.0, + "urban_runoff": 0.0, + } + ) + state = pd.DataFrame(data={"node_id": [2], "level": 0.5}) + basin = ribasim.Basin(profile=profile, static=static, state=state) + + # Setup flow boundary + flow_boundary = ribasim.FlowBoundary( + static=pd.DataFrame(data={"node_id": [1], "flow_rate": 1.0}) + ) + + # Setup allocation level control + allocation_level_control = ribasim.AllocationLevelControl( + static=pd.DataFrame(data={"node_id": [4], "priority": 1, "target_level": 1.0}) + ) + + # Setup user + user = ribasim.User( + static=pd.DataFrame( + data={ + "node_id": [3], + "priority": [2], + "demand": [1.0], + "return_factor": [0.0], + "min_level": [0.2], + } + ) + ) + + # Setup allocation + allocation = ribasim.Allocation(use_allocation=True, timestep=86400) + + model = ribasim.Model( + network=ribasim.Network(node=node, edge=edge), + basin=basin, + flow_boundary=flow_boundary, + allocation_level_control=allocation_level_control, + user=user, + allocation=allocation, + starttime="2020-01-01 00:00:00", + endtime="2020-03-01 00:00:00", + ) + + return model From 9812528efe0683fcd59c9bf693b2ab81bd4293b5 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 6 Feb 2024 16:13:29 +0100 Subject: [PATCH 06/37] Add basin flow variables and include in flow conservation constraints --- core/src/allocation.jl | 39 ++++++++++++++++++++++++++++++++++++--- docs/core/allocation.qmd | 6 ++---- 2 files changed, 38 insertions(+), 7 deletions(-) diff --git a/core/src/allocation.jl b/core/src/allocation.jl index 2c4b8e2dd..9aad3a634 100644 --- a/core/src/allocation.jl +++ b/core/src/allocation.jl @@ -401,6 +401,24 @@ function add_variables_flow!( return nothing end +""" +Add the variables for supply/demand of a basin to the problem. +The variable indices are the node_ids of the basins in the subnetwork. +""" +function add_variables_basin!( + problem::JuMP.Model, + p::Parameters, + allocation_network_id::Int, +)::Nothing + (; graph) = p + node_ids_basin = [ + node_id for node_id in graph[].node_ids[allocation_network_id] if + graph[node_id].type == :basin + ] + problem[:F_basin] = JuMP.@variable(problem, F_basin[node_id = node_ids_basin,]) + 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, @@ -523,6 +541,18 @@ function outflow_ids_allocation(graph::MetaGraph, node_id::NodeID) return outflow_ids end +function get_basin_flow( + problem::JuMP.Model, + node_id::NodeID, +)::Union{JuMP.VariableRef, Float64} + F_basin = problem[:F_basin] + return if node_id in only(F_basin.axes) + F_basin[node_id] + else + 0.0 + end +end + """ Add the flow conservation constraints to the allocation problem. The constraint indices are user node IDs. @@ -535,8 +565,9 @@ function add_constraints_flow_conservation!( p::Parameters, allocation_network_id::Int, )::Nothing - (; graph, allocation) = p + (; graph) = p F = problem[:F] + F_basin = problem[:F_basin] node_ids = graph[].node_ids[allocation_network_id] node_ids_conservation = [node_id for node_id in node_ids if graph[node_id].type == :basin] @@ -551,10 +582,11 @@ 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) - ]), + ]) + get_basin_flow(problem, node_id), base_name = "flow_conservation", ) return nothing @@ -732,6 +764,7 @@ 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, config) # TODO: Add variables for allocation to basins diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index e09d11765..074ff6b86 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -107,10 +107,8 @@ The optimization problem for a subnetwork is a linear optimization problem consi There are 2 types of variable whose value has to be determined to solve the allocation problem: - The flows $F \in \mathbb{R}_{\ge 0}^{n'\times n'}$ over the edges in the allocation graph; -- The allocations to the basins -$$ - A^\text{basin}_{i} \ge 0, \quad B_S. -$$ +- The flows $A^\text{basin}_{i}$ for all $ \quad B_S$ supplied (positive) or consumed (negative) by the basins. + :::{.callout-note} Currently the basin allocations are not taken into account in the implementation. From c24143f0012d4702df33c80dccef44f81b99cea2 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 7 Feb 2024 09:24:57 +0100 Subject: [PATCH 07/37] Accepted level window instead of single level --- core/src/create.jl | 11 +++++++---- core/src/solve.jl | 6 ++++-- core/src/validation.jl | 6 ++++-- docs/core/allocation.qmd | 4 ++-- docs/schema/AllocationLevelControlStatic.schema.json | 9 +++++++-- docs/schema/AllocationLevelControlTime.schema.json | 9 +++++++-- python/ribasim/ribasim/models.py | 6 ++++-- .../ribasim_testmodels/allocation.py | 4 +++- 8 files changed, 38 insertions(+), 17 deletions(-) diff --git a/core/src/create.jl b/core/src/create.jl index 4858856b7..1f75f0339 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -768,19 +768,22 @@ function AllocationLevelControl(db::DB, config::Config)::AllocationLevelControl "AllocationLevelControl"; static, time, - time_interpolatables = [:target_level], + time_interpolatables = [:min_level, :max_level], ) if !valid error("Errors occurred when parsing AllocationLevelControl data.") end - target_level = zeros(length(parsed_parameters.node_id)) + min_level = zeros(length(parsed_parameters.node_id)) + max_level = zeros(length(parsed_parameters.node_id)) return AllocationLevelControl( NodeID.(parsed_parameters.node_id), - target_level, - parsed_parameters.target_level, + min_level, + max_level, + parsed_parameters.min_level, + parsed_parameters.max_level, parsed_parameters.priority, ) end diff --git a/core/src/solve.jl b/core/src/solve.jl index eec444f4a..8a3f70429 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -476,8 +476,10 @@ end struct AllocationLevelControl node_id::Vector{NodeID} - target_level::Vector{Float64} - target_level_itp::Vector{LinearInterpolation} + min_level::Vector{Float64} + max_level::Vector{Float64} + min_level_itp::Vector{LinearInterpolation} + max_level_itp::Vector{LinearInterpolation} priority::Vector{Int} end diff --git a/core/src/validation.jl b/core/src/validation.jl index b3cbf26ab..47be4ebb1 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -338,14 +338,16 @@ end @version AllocationLevelControlStaticV1 begin node_id::Int - target_level::Float64 + min_level::Float64 + max_level::Float64 priority::Int end @version AllocationLevelControlTimeV1 begin node_id::Int time::DateTime - target_level::Float64 + min_level::Float64 + max_level::Float64 priority::Int end diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index 074ff6b86..44fe0afae 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -107,7 +107,7 @@ The optimization problem for a subnetwork is a linear optimization problem consi There are 2 types of variable whose value has to be determined to solve the allocation problem: - The flows $F \in \mathbb{R}_{\ge 0}^{n'\times n'}$ over the edges in the allocation graph; -- The flows $A^\text{basin}_{i}$ for all $ \quad B_S$ supplied (positive) or consumed (negative) by the basins. +- The flows $A^\text{basin}_{i}$ for all $i \in B_S$ supplied (positive) or consumed (negative) by the basins. :::{.callout-note} @@ -159,7 +159,7 @@ 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} \le A^\text{basin}_k + \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: diff --git a/docs/schema/AllocationLevelControlStatic.schema.json b/docs/schema/AllocationLevelControlStatic.schema.json index b04fc69a7..18dcfa84c 100644 --- a/docs/schema/AllocationLevelControlStatic.schema.json +++ b/docs/schema/AllocationLevelControlStatic.schema.json @@ -9,7 +9,11 @@ "format": "default", "type": "integer" }, - "target_level": { + "min_level": { + "format": "double", + "type": "number" + }, + "max_level": { "format": "double", "type": "number" }, @@ -26,7 +30,8 @@ }, "required": [ "node_id", - "target_level", + "min_level", + "max_level", "priority" ] } diff --git a/docs/schema/AllocationLevelControlTime.schema.json b/docs/schema/AllocationLevelControlTime.schema.json index 5478922de..962683800 100644 --- a/docs/schema/AllocationLevelControlTime.schema.json +++ b/docs/schema/AllocationLevelControlTime.schema.json @@ -13,7 +13,11 @@ "format": "date-time", "type": "string" }, - "target_level": { + "min_level": { + "format": "double", + "type": "number" + }, + "max_level": { "format": "double", "type": "number" }, @@ -31,7 +35,8 @@ "required": [ "node_id", "time", - "target_level", + "min_level", + "max_level", "priority" ] } diff --git a/python/ribasim/ribasim/models.py b/python/ribasim/ribasim/models.py index 3e7faa7bb..2227774ff 100644 --- a/python/ribasim/ribasim/models.py +++ b/python/ribasim/ribasim/models.py @@ -12,7 +12,8 @@ class AllocationLevelControlStatic(BaseModel): node_id: int - target_level: float + min_level: float + max_level: float priority: int remarks: str = Field("", description="a hack for pandera") @@ -20,7 +21,8 @@ class AllocationLevelControlStatic(BaseModel): class AllocationLevelControlTime(BaseModel): node_id: int time: datetime - target_level: float + min_level: float + max_level: float priority: int remarks: str = Field("", description="a hack for pandera") diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 38dc458b5..580b2a5bb 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -1808,7 +1808,9 @@ def allocation_level_control_model(): # Setup allocation level control allocation_level_control = ribasim.AllocationLevelControl( - static=pd.DataFrame(data={"node_id": [4], "priority": 1, "target_level": 1.0}) + static=pd.DataFrame( + data={"node_id": [4], "priority": 1, "min_level": 1.0, "max_level": 1.0} + ) ) # Setup user From c714de93a72b1c31dc68a9d5864d3f91634ecc8e Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 7 Feb 2024 11:30:40 +0100 Subject: [PATCH 08/37] Towards basin flow constraints --- core/src/allocation.jl | 68 +++++++++++++++++++++++++++++++++++----- core/src/bmi.jl | 6 ++-- docs/core/allocation.qmd | 6 ++-- 3 files changed, 68 insertions(+), 12 deletions(-) diff --git a/core/src/allocation.jl b/core/src/allocation.jl index 9aad3a634..fd7702708 100644 --- a/core/src/allocation.jl +++ b/core/src/allocation.jl @@ -750,6 +750,25 @@ function add_constraints_fractional_flow!( return nothing end +""" +Add the basin flow constraints to the allocation problem. +The constraint indices are the basin node IDs. +NOTE: negative basin flow means flow out of the basin. + +Constraint: +flow basin >= - basin capacity +""" +function add_constraints_basin_flow!(problem::JuMP.Model)::Nothing + F_basin = problem[:F_basin] + problem[:basin_flow] = JuMP.@constraint( + problem, + [node_id = only(F_basin.axes)], + F_basin[node_id] >= 0.0, + base_name = "basin_flow" + ) + return nothing +end + """ Construct the allocation problem for the current subnetwork as a JuMP.jl model. """ @@ -766,7 +785,6 @@ function allocation_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, config) - # TODO: Add variables for allocation to basins # Add constraints to problem add_constraints_capacity!(problem, capacity, p, allocation_network_id) @@ -775,6 +793,7 @@ function allocation_problem( add_constraints_user_returnflow!(problem, p, allocation_network_id) add_constraints_absolute_value!(problem, p, allocation_network_id, config) add_constraints_fractional_flow!(problem, p, allocation_network_id) + add_constraints_basin_flow!(problem) return problem end @@ -1105,6 +1124,43 @@ function adjust_edge_capacities!( end end +function get_basin_capacity( + allocation_model::AllocationModel, + u::ComponentVector, + p::Parameters, + node_id::NodeID, +)::Float64 + (; graph, user, basin) = p + (; Δt_allocation) = allocation_model + influx = get_flow(graph, node_id, 0.0) + basin_idx = findsorted(user.node_id, node_id) + storage_basin = u.storage[basin_idx] + level_max = ... + storage_max = get_storage_from_level(basin, basin_idx, level_max) + + return min(0.0, ()) +end + +function adjust_basin_capacities!( + allocation_model::AllocationModel, + u::ComponentVector, + p::Parameters, + priority_idx::Int, +)::Nothing + (; problem) = allocation_model + constraints_basin_flow = problem[:basin_flow] + + if priority_idx == 1 + for node_id in only(constraints_basin_flow.axes) + JuMP.set_normalized_rhs( + constraints_basin_flow[node_id], + -get_basin_capacity(allocation_model, u), + ) + end + end + return nothing +end + """ Save the allocation flows per physical edge. """ @@ -1146,7 +1202,8 @@ and flows, solve the allocation problem and assign the results to the users. function allocate!( p::Parameters, allocation_model::AllocationModel, - t::Float64; + t::Float64, + u::ComponentVector; collect_demands::Bool = false, )::Nothing (; user, allocation) = p @@ -1154,11 +1211,6 @@ function allocate!( (; priorities) = user (; subnetwork_demands) = allocation - # TODO: Compute basin flow from vertical fluxes and basin volume. - # Set as basin demand if the net flow is negative, set as source - # in the flow_conservation constraints if the net flow is positive. - # Solve this as a separate problem before the priorities below - main_network_source_edges = get_main_network_connections(p, allocation_network_id) if collect_demands @@ -1176,6 +1228,8 @@ function allocate!( # or set edge capacities if priority_idx = 1 adjust_edge_capacities!(allocation_model, p, priority_idx) + adjust_basin_capacities!(allocation_model, u, p, priority_idx) + # Set the objective depending on the demands # A new objective function is set instead of modifying the coefficients # of an existing objective function because this is not supported for diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 8b6df4d77..f3d8c5889 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -607,14 +607,14 @@ end "Solve the allocation problem for all users and assign allocated abstractions to user nodes." function update_allocation!(integrator)::Nothing - (; p, t) = integrator + (; p, t, u) = integrator (; allocation) = p (; allocation_models) = allocation # If a main network is present, collect demands of subnetworks if has_main_network(allocation) for allocation_model in Iterators.drop(allocation_models, 1) - allocate!(p, allocation_model, t; collect_demands = true) + allocate!(p, allocation_model, t, u; collect_demands = true) end end @@ -622,7 +622,7 @@ function update_allocation!(integrator)::Nothing # If a main network is present this is solved first, # which provides allocation to the subnetworks for allocation_model in allocation_models - allocate!(p, allocation_model, t) + allocate!(p, allocation_model, t, u) end end diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index 44fe0afae..afed1f0f5 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -54,6 +54,8 @@ $$ \phi_i(t), \quad \forall i \in B_S. $$ +We consider fluxes into the basin to be positive and out of the basin to be negative. + Secondly, there is the available water in each basin above the minimum level $l_{\min,i}$ corresponding to a minimum storage $s_{\min,i}$: $$ u_i(t)-s_{\min,i}, \quad \forall i \in B_S. @@ -107,7 +109,7 @@ The optimization problem for a subnetwork is a linear optimization problem consi There are 2 types of variable whose value has to be determined to solve the allocation problem: - The flows $F \in \mathbb{R}_{\ge 0}^{n'\times n'}$ over the edges in the allocation graph; -- The flows $A^\text{basin}_{i}$ for all $i \in B_S$ supplied (positive) or consumed (negative) by the basins. +- The flows $F^\text{basin}_{i}$ for all $i \in B_S$ supplied (positive) or consumed (negative) by the basins. :::{.callout-note} @@ -159,7 +161,7 @@ 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 A^\text{basin}_k + \sum_{i=1}^{n'} F_{ik}, \quad \forall k \in B_S . + \sum_{j=1}^{n'} F_{kj} \le F^\text{basin}_k + \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: From 8a89906f6e565eacf9d6cbbf2ee1df4826522806 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 7 Feb 2024 13:04:38 +0100 Subject: [PATCH 09/37] Complete basin capacity constraints --- core/src/allocation.jl | 32 +++++++++++++++++++++++--------- core/src/create.jl | 5 ----- core/src/solve.jl | 8 +++----- 3 files changed, 26 insertions(+), 19 deletions(-) diff --git a/core/src/allocation.jl b/core/src/allocation.jl index fd7702708..fecc56d52 100644 --- a/core/src/allocation.jl +++ b/core/src/allocation.jl @@ -1128,36 +1128,50 @@ function get_basin_capacity( allocation_model::AllocationModel, u::ComponentVector, p::Parameters, + t::Float64, node_id::NodeID, )::Float64 - (; graph, user, basin) = p + (; graph, basin, allocation_level_control) = p (; Δt_allocation) = allocation_model influx = get_flow(graph, node_id, 0.0) - basin_idx = findsorted(user.node_id, node_id) + _, basin_idx = id_index(basin.node_id, node_id) storage_basin = u.storage[basin_idx] - level_max = ... + allocation_level_control_node_id = + first(inneighbor_labels_type(graph, node_id, EdgeType.control)) + allocation_level_control_idx = + findsorted(allocation_level_control.node_id, allocation_level_control_node_id) + level_max = allocation_level_control.max_level[allocation_level_control_idx](t) storage_max = get_storage_from_level(basin, basin_idx, level_max) - return min(0.0, ()) + return max(0.0, (storage_basin - storage_max) / Δt_allocation + influx) end function adjust_basin_capacities!( allocation_model::AllocationModel, u::ComponentVector, p::Parameters, + t::Float64, priority_idx::Int, )::Nothing (; problem) = allocation_model constraints_basin_flow = problem[:basin_flow] + F_basin = problem[:F_basin] - if priority_idx == 1 - for node_id in only(constraints_basin_flow.axes) + for node_id in only(constraints_basin_flow.axes) + constraint = constraints_basin_flow[node_id] + if priority_idx == 1 JuMP.set_normalized_rhs( - constraints_basin_flow[node_id], - -get_basin_capacity(allocation_model, u), + constraint, + -get_basin_capacity(allocation_model, u, p, t, node_id), + ) + else + JuMP.set_normalized_rhs( + constraint, + JuMP.normalized_rhs(constraint) - JuMP.value(F_basin[node_id]), ) end end + return nothing end @@ -1228,7 +1242,7 @@ function allocate!( # or set edge capacities if priority_idx = 1 adjust_edge_capacities!(allocation_model, p, priority_idx) - adjust_basin_capacities!(allocation_model, u, p, priority_idx) + adjust_basin_capacities!(allocation_model, u, p, t, priority_idx) # Set the objective depending on the demands # A new objective function is set instead of modifying the coefficients diff --git a/core/src/create.jl b/core/src/create.jl index 1f75f0339..5c39a83fd 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -775,13 +775,8 @@ function AllocationLevelControl(db::DB, config::Config)::AllocationLevelControl error("Errors occurred when parsing AllocationLevelControl data.") end - min_level = zeros(length(parsed_parameters.node_id)) - max_level = zeros(length(parsed_parameters.node_id)) - return AllocationLevelControl( NodeID.(parsed_parameters.node_id), - min_level, - max_level, parsed_parameters.min_level, parsed_parameters.max_level, parsed_parameters.priority, diff --git a/core/src/solve.jl b/core/src/solve.jl index 8a3f70429..389185494 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -476,10 +476,8 @@ end struct AllocationLevelControl node_id::Vector{NodeID} - min_level::Vector{Float64} - max_level::Vector{Float64} - min_level_itp::Vector{LinearInterpolation} - max_level_itp::Vector{LinearInterpolation} + min_level::Vector{LinearInterpolation} + max_level::Vector{LinearInterpolation} priority::Vector{Int} end @@ -1174,7 +1172,7 @@ function formulate_du!( basin::Basin, storage::AbstractVector, )::Nothing - (; flow_vertical_dict, flow_vertical) = graph[] + (; flow_vertical) = graph[] flow_vertical = get_tmp(flow_vertical, storage) # loop over basins # subtract all outgoing flows From 25cba1ca8639bdddb1685fdf5da18146dea7315a Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 7 Feb 2024 13:30:03 +0100 Subject: [PATCH 10/37] Update testmodel --- .../ribasim_testmodels/allocation.py | 33 ++++++++++--------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 580b2a5bb..1f3180bc0 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -1747,18 +1747,19 @@ def allocation_level_control_model(): (1.0, 0.0), # 2: Basin (2.0, 0.0), # 3: User (1.0, -1.0), # 4: AllocationLevelControl + (2.0, -1.0), # 5: Basin ] ) node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1]) - node_type = ["FlowBoundary", "Basin", "User", "AllocationLevelControl"] + node_type = ["FlowBoundary", "Basin", "User", "AllocationLevelControl", "Basin"] # Make sure the feature id starts at 1: explicitly give an index. node = ribasim.Node( df=gpd.GeoDataFrame( data={ "type": node_type, - "allocation_network_id": 4 * [1], + "allocation_network_id": 5 * [2], }, index=pd.Index(np.arange(len(xy)) + 1, name="fid"), geometry=node_xy, @@ -1767,10 +1768,10 @@ def allocation_level_control_model(): ) # Setup the edges: - from_id = np.array([1, 2, 4, 3]) - to_id = np.array([2, 3, 2, 2]) - edge_type = ["flow", "flow", "control", "flow"] - allocation_network_id = [1, None, None, None] + from_id = np.array([1, 2, 4, 3, 4]) + to_id = np.array([2, 3, 2, 5, 5]) + edge_type = ["flow", "flow", "control", "flow", "control"] + allocation_network_id = [1, None, None, None, None] lines = node.geometry_from_connectivity(from_id.tolist(), to_id.tolist()) edge = ribasim.Edge( @@ -1787,23 +1788,25 @@ def allocation_level_control_model(): ) # Setup basin - profile = pd.DataFrame(data={"node_id": 2, "area": 1e3, "level": [0.0, 1.0]}) + profile = pd.DataFrame( + data={"node_id": [2, 2, 5, 5], "area": 1e3, "level": [0.0, 1.0, 0.0, 1.0]} + ) static = pd.DataFrame( data={ - "node_id": [2], - "drainage": 0.0, - "potential_evaporation": 0.0, - "infiltration": 0.0, - "precipitation": 0.0, + "node_id": [2, 5], + "drainage": 0.1e-3, + "potential_evaporation": 0.2e-3, + "infiltration": 0.3e-3, + "precipitation": 0.4e-3, "urban_runoff": 0.0, } ) - state = pd.DataFrame(data={"node_id": [2], "level": 0.5}) + state = pd.DataFrame(data={"node_id": [2, 5], "level": 0.5}) basin = ribasim.Basin(profile=profile, static=static, state=state) # Setup flow boundary flow_boundary = ribasim.FlowBoundary( - static=pd.DataFrame(data={"node_id": [1], "flow_rate": 1.0}) + static=pd.DataFrame(data={"node_id": [1], "flow_rate": 1e-3}) ) # Setup allocation level control @@ -1820,7 +1823,7 @@ def allocation_level_control_model(): "node_id": [3], "priority": [2], "demand": [1.0], - "return_factor": [0.0], + "return_factor": [0.2], "min_level": [0.2], } ) From 821c57486d8e78fd1c1e0ca91dde5680284ee5d9 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 7 Feb 2024 13:41:47 +0100 Subject: [PATCH 11/37] Split allocation code into inititialization and optimization --- core/src/Ribasim.jl | 3 +- .../src/{allocation.jl => allocation_init.jl} | 444 ------------------ core/src/allocation_optim.jl | 443 +++++++++++++++++ 3 files changed, 445 insertions(+), 445 deletions(-) rename core/src/{allocation.jl => allocation_init.jl} (65%) create mode 100644 core/src/allocation_optim.jl diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index ca675a366..819db25a3 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -80,7 +80,8 @@ include("solve.jl") include("config.jl") using .config include("logging.jl") -include("allocation.jl") +include("allocation_init.jl") +include("allocation_optim.jl") include("utils.jl") include("lib.jl") include("io.jl") diff --git a/core/src/allocation.jl b/core/src/allocation_init.jl similarity index 65% rename from core/src/allocation.jl rename to core/src/allocation_init.jl index fecc56d52..77a923633 100644 --- a/core/src/allocation.jl +++ b/core/src/allocation_init.jl @@ -831,447 +831,3 @@ function AllocationModel( Δt_allocation, ) end - -""" -Add a term to the expression of the objective function corresponding to -the demand of a user. -""" -function add_user_term!( - ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, - edge::Tuple{NodeID, NodeID}, - objective_type::Symbol, - demand::Float64, - model::AllocationModel, -)::Nothing - (; problem) = model - F = problem[:F] - F_edge = F[edge] - node_id_user = edge[2] - - if objective_type == :quadratic_absolute - # Objective function ∑ (F - d)^2 - JuMP.add_to_expression!(ex, 1, F_edge, F_edge) - JuMP.add_to_expression!(ex, -2 * demand, F_edge) - JuMP.add_to_expression!(ex, demand^2) - - elseif objective_type == :quadratic_relative - # Objective function ∑ (1 - F/d)^2 - if demand ≈ 0 - return nothing - end - JuMP.add_to_expression!(ex, 1.0 / demand^2, F_edge, F_edge) - JuMP.add_to_expression!(ex, -2.0 / demand, F_edge) - JuMP.add_to_expression!(ex, 1.0) - - elseif objective_type == :linear_absolute - # Objective function ∑ |F - d| - JuMP.set_normalized_rhs(problem[:abs_positive][node_id_user], -demand) - JuMP.set_normalized_rhs(problem[:abs_negative][node_id_user], demand) - - elseif objective_type == :linear_relative - # Objective function ∑ |1 - F/d| - JuMP.set_normalized_coefficient( - problem[:abs_positive][node_id_user], - F_edge, - iszero(demand) ? 0 : 1 / demand, - ) - JuMP.set_normalized_coefficient( - problem[:abs_negative][node_id_user], - F_edge, - iszero(demand) ? 0 : -1 / demand, - ) - else - error("Invalid allocation objective type $objective_type.") - end - return nothing -end - -""" -Set the objective for the given priority. -For an objective with absolute values this also involves adjusting constraints. -""" -function set_objective_priority!( - allocation_model::AllocationModel, - p::Parameters, - t::Float64, - priority_idx::Int, -)::Nothing - (; objective_type, problem, allocation_network_id) = allocation_model - (; graph, user, allocation) = p - (; demand, node_id) = user - (; main_network_connections, subnetwork_demands) = allocation - edge_ids = graph[].edge_ids[allocation_network_id] - - F = problem[:F] - if objective_type in [:quadratic_absolute, :quadratic_relative] - ex = JuMP.QuadExpr() - elseif objective_type in [:linear_absolute, :linear_relative] - ex = JuMP.AffExpr() - 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 - end - - # Terms for user nodes - for edge_id in edge_ids - node_id_user = edge_id[2] - if graph[node_id_user].type != :user - continue - end - - user_idx = findsorted(node_id, node_id_user) - d = demand[user_idx][priority_idx](t) - 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 -end - -""" -Assign the allocations to the users as determined by the solution of the allocation problem. -""" -function assign_allocations!( - allocation_model::AllocationModel, - p::Parameters, - t::Float64, - priority_idx::Int; - collect_demands::Bool = false, -)::Nothing - (; problem, allocation_network_id) = allocation_model - (; graph, user, allocation) = p - (; - subnetwork_demands, - subnetwork_allocateds, - allocation_network_ids, - main_network_connections, - ) = allocation - (; record) = user - edge_ids = graph[].edge_ids[allocation_network_id] - main_network_source_edges = get_main_network_connections(p, allocation_network_id) - F = problem[:F] - for edge_id in edge_ids - # If this edge is a source edge from the main network to a subnetwork, - # and demands are being collected, add its flow to the demand of this edge - if collect_demands && - graph[edge_id...].allocation_network_id_source == allocation_network_id && - edge_id ∈ main_network_source_edges - allocated = JuMP.value(F[edge_id]) - subnetwork_demands[edge_id][priority_idx] += allocated - end - - user_node_id = edge_id[2] - - if graph[user_node_id].type == :user - allocated = JuMP.value(F[edge_id]) - user_idx = findsorted(user.node_id, user_node_id) - user.allocated[user_idx][priority_idx] = allocated - - # Save allocations to record - push!(record.time, t) - push!(record.allocation_network_id, allocation_model.allocation_network_id) - push!(record.user_node_id, Int(user_node_id)) - push!(record.priority, user.priorities[priority_idx]) - push!(record.demand, user.demand[user_idx][priority_idx](t)) - push!(record.allocated, allocated) - # TODO: This is now the last abstraction before the allocation update, - # should be the average abstraction since the last allocation solve - push!( - record.abstracted, - get_flow(graph, inflow_id(graph, user_node_id), user_node_id, 0), - ) - end - end - - # Write the flows to the subnetworks as allocated flows - # in the allocation object - if is_main_network(allocation_network_id) - for (allocation_network_id, main_network_source_edges) in - zip(allocation_network_ids, main_network_connections) - if is_main_network(allocation_network_id) - continue - end - for edge_id in main_network_source_edges - subnetwork_allocateds[edge_id][priority_idx] = JuMP.value(F[edge_id]) - end - end - end - return nothing -end - -""" -Adjust the source flows. -""" -function adjust_source_capacities!( - allocation_model::AllocationModel, - p::Parameters, - priority_idx::Int; - collect_demands::Bool = false, -)::Nothing - (; problem) = allocation_model - (; graph, allocation) = p - (; allocation_network_id) = allocation_model - (; subnetwork_allocateds) = allocation - edge_ids = graph[].edge_ids[allocation_network_id] - source_constraints = problem[:source] - F = problem[:F] - - main_network_source_edges = get_main_network_connections(p, allocation_network_id) - - 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 priority_idx == 1 - # If the optimization was just started, i.e. sources have to be reset - if edge_id in main_network_source_edges - if collect_demands - # Set the source capacity to effectively unlimited if subnetwork demands are being collected - source_capacity = Inf - else - # Set the source capacity to the value allocated to the subnetwork over this edge - source_capacity = subnetwork_allocateds[edge_id][priority_idx] - end - else - # Reset the source to the current flow from the physical layer. - source_capacity = get_flow(graph, edge_id..., 0) - end - JuMP.set_normalized_rhs( - source_constraints[edge_id], - # It is assumed that the allocation procedure does not have to be differentiated. - source_capacity, - ) - else - # Subtract the allocated flow from the source. - JuMP.set_normalized_rhs( - source_constraints[edge_id], - JuMP.normalized_rhs(source_constraints[edge_id]) - - JuMP.value(F[edge_id]), - ) - end - end - end - return nothing -end - -""" -Set the values of the edge capacities. 2 cases: -- Before the first allocation solve, set the edge capacities to their full capacity; -- Before an allocation solve, subtract the flow used by allocation for the previous priority - from the edge capacities. -""" -function adjust_edge_capacities!( - allocation_model::AllocationModel, - p::Parameters, - priority_idx::Int, -)::Nothing - (; graph) = p - (; problem, capacity, allocation_network_id) = allocation_model - edge_ids = graph[].edge_ids[allocation_network_id] - constraints_capacity = problem[:capacity] - F = problem[:F] - - main_network_source_edges = get_main_network_connections(p, allocation_network_id) - - for edge_id in edge_ids - c = capacity[edge_id...] - - # These edges have no capacity constraints: - # - With infinite capacity - # - Being a source from the main network to a subnetwork - if isinf(c) || edge_id ∈ main_network_source_edges - continue - end - - if priority_idx == 1 - # Before the first allocation solve, set the edge capacities to their full capacity - JuMP.set_normalized_rhs(constraints_capacity[edge_id], c) - else - # Before an allocation solve, subtract the flow used by allocation for the previous priority - # from the edge capacities - JuMP.set_normalized_rhs( - constraints_capacity[edge_id], - JuMP.normalized_rhs(constraints_capacity[edge_id]) - JuMP.value(F[edge_id]), - ) - end - end -end - -function get_basin_capacity( - allocation_model::AllocationModel, - u::ComponentVector, - p::Parameters, - t::Float64, - node_id::NodeID, -)::Float64 - (; graph, basin, allocation_level_control) = p - (; Δt_allocation) = allocation_model - influx = get_flow(graph, node_id, 0.0) - _, basin_idx = id_index(basin.node_id, node_id) - storage_basin = u.storage[basin_idx] - allocation_level_control_node_id = - first(inneighbor_labels_type(graph, node_id, EdgeType.control)) - allocation_level_control_idx = - findsorted(allocation_level_control.node_id, allocation_level_control_node_id) - level_max = allocation_level_control.max_level[allocation_level_control_idx](t) - storage_max = get_storage_from_level(basin, basin_idx, level_max) - - return max(0.0, (storage_basin - storage_max) / Δt_allocation + influx) -end - -function adjust_basin_capacities!( - allocation_model::AllocationModel, - u::ComponentVector, - p::Parameters, - t::Float64, - priority_idx::Int, -)::Nothing - (; problem) = allocation_model - constraints_basin_flow = problem[:basin_flow] - F_basin = problem[:F_basin] - - for node_id in only(constraints_basin_flow.axes) - constraint = constraints_basin_flow[node_id] - if priority_idx == 1 - JuMP.set_normalized_rhs( - constraint, - -get_basin_capacity(allocation_model, u, p, t, node_id), - ) - else - JuMP.set_normalized_rhs( - constraint, - JuMP.normalized_rhs(constraint) - JuMP.value(F_basin[node_id]), - ) - end - end - - return nothing -end - -""" -Save the allocation flows per physical edge. -""" -function save_allocation_flows!( - p::Parameters, - t::Float64, - allocation_model::AllocationModel, - priority::Int, - collect_demands::Bool, -)::Nothing - (; problem, allocation_network_id) = allocation_model - (; allocation, graph) = p - (; record) = allocation - F = problem[:F] - - for allocation_edge in first(F.axes) - flow = JuMP.value(F[allocation_edge]) - edge_metadata = graph[allocation_edge...] - (; node_ids) = edge_metadata - - for i in eachindex(node_ids)[1:(end - 1)] - push!(record.time, t) - push!(record.edge_id, edge_metadata.id) - push!(record.from_node_id, node_ids[i]) - push!(record.to_node_id, node_ids[i + 1]) - push!(record.allocation_network_id, allocation_network_id) - push!(record.priority, priority) - push!(record.flow, flow) - push!(record.collect_demands, collect_demands) - end - end - return nothing -end - -""" -Update the allocation optimization problem for the given subnetwork with the problem state -and flows, solve the allocation problem and assign the results to the users. -""" -function allocate!( - p::Parameters, - allocation_model::AllocationModel, - t::Float64, - u::ComponentVector; - collect_demands::Bool = false, -)::Nothing - (; user, allocation) = p - (; problem, allocation_network_id) = allocation_model - (; priorities) = user - (; subnetwork_demands) = allocation - - main_network_source_edges = get_main_network_connections(p, allocation_network_id) - - if collect_demands - for main_network_connection in keys(subnetwork_demands) - if main_network_connection in main_network_source_edges - subnetwork_demands[main_network_connection] .= 0.0 - end - end - end - - for priority_idx in eachindex(priorities) - adjust_source_capacities!(allocation_model, p, priority_idx; collect_demands) - - # Subtract the flows used by the allocation of the previous priority from the capacities of the edges - # or set edge capacities if priority_idx = 1 - adjust_edge_capacities!(allocation_model, p, priority_idx) - - adjust_basin_capacities!(allocation_model, u, p, t, priority_idx) - - # Set the objective depending on the demands - # A new objective function is set instead of modifying the coefficients - # of an existing objective function because this is not supported for - # quadratic terms: - # https://jump.dev/JuMP.jl/v1.16/manual/objective/#Modify-an-objective-coefficient - set_objective_priority!(allocation_model, p, t, priority_idx) - - # Solve the allocation problem for this priority - JuMP.optimize!(problem) - @debug JuMP.solution_summary(problem) - if JuMP.termination_status(problem) !== JuMP.OPTIMAL - (; allocation_network_id) = allocation_model - priority = priorities[priority_index] - error( - "Allocation of subnetwork $allocation_network_id, priority $priority coudn't find optimal solution.", - ) - end - - # Assign the allocations to the users for this priority - assign_allocations!(allocation_model, p, t, priority_idx; collect_demands) - - # Save the flows over all edges in the subnetwork - save_allocation_flows!( - p, - t, - allocation_model, - priorities[priority_idx], - collect_demands, - ) - end -end diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl new file mode 100644 index 000000000..ca218b76e --- /dev/null +++ b/core/src/allocation_optim.jl @@ -0,0 +1,443 @@ +""" +Add a term to the expression of the objective function corresponding to +the demand of a user. +""" +function add_user_term!( + ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, + edge::Tuple{NodeID, NodeID}, + objective_type::Symbol, + demand::Float64, + model::AllocationModel, +)::Nothing + (; problem) = model + F = problem[:F] + F_edge = F[edge] + node_id_user = edge[2] + + if objective_type == :quadratic_absolute + # Objective function ∑ (F - d)^2 + JuMP.add_to_expression!(ex, 1, F_edge, F_edge) + JuMP.add_to_expression!(ex, -2 * demand, F_edge) + JuMP.add_to_expression!(ex, demand^2) + + elseif objective_type == :quadratic_relative + # Objective function ∑ (1 - F/d)^2 + if demand ≈ 0 + return nothing + end + JuMP.add_to_expression!(ex, 1.0 / demand^2, F_edge, F_edge) + JuMP.add_to_expression!(ex, -2.0 / demand, F_edge) + JuMP.add_to_expression!(ex, 1.0) + + elseif objective_type == :linear_absolute + # Objective function ∑ |F - d| + JuMP.set_normalized_rhs(problem[:abs_positive][node_id_user], -demand) + JuMP.set_normalized_rhs(problem[:abs_negative][node_id_user], demand) + + elseif objective_type == :linear_relative + # Objective function ∑ |1 - F/d| + JuMP.set_normalized_coefficient( + problem[:abs_positive][node_id_user], + F_edge, + iszero(demand) ? 0 : 1 / demand, + ) + JuMP.set_normalized_coefficient( + problem[:abs_negative][node_id_user], + F_edge, + iszero(demand) ? 0 : -1 / demand, + ) + else + error("Invalid allocation objective type $objective_type.") + end + return nothing +end + +""" +Set the objective for the given priority. +For an objective with absolute values this also involves adjusting constraints. +""" +function set_objective_priority!( + allocation_model::AllocationModel, + p::Parameters, + t::Float64, + priority_idx::Int, +)::Nothing + (; objective_type, problem, allocation_network_id) = allocation_model + (; graph, user, allocation) = p + (; demand, node_id) = user + (; main_network_connections, subnetwork_demands) = allocation + edge_ids = graph[].edge_ids[allocation_network_id] + + F = problem[:F] + if objective_type in [:quadratic_absolute, :quadratic_relative] + ex = JuMP.QuadExpr() + elseif objective_type in [:linear_absolute, :linear_relative] + ex = JuMP.AffExpr() + 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 + end + + # Terms for user nodes + for edge_id in edge_ids + node_id_user = edge_id[2] + if graph[node_id_user].type != :user + continue + end + + user_idx = findsorted(node_id, node_id_user) + d = demand[user_idx][priority_idx](t) + 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 +end + +""" +Assign the allocations to the users as determined by the solution of the allocation problem. +""" +function assign_allocations!( + allocation_model::AllocationModel, + p::Parameters, + t::Float64, + priority_idx::Int; + collect_demands::Bool = false, +)::Nothing + (; problem, allocation_network_id) = allocation_model + (; graph, user, allocation) = p + (; + subnetwork_demands, + subnetwork_allocateds, + allocation_network_ids, + main_network_connections, + ) = allocation + (; record) = user + edge_ids = graph[].edge_ids[allocation_network_id] + main_network_source_edges = get_main_network_connections(p, allocation_network_id) + F = problem[:F] + for edge_id in edge_ids + # If this edge is a source edge from the main network to a subnetwork, + # and demands are being collected, add its flow to the demand of this edge + if collect_demands && + graph[edge_id...].allocation_network_id_source == allocation_network_id && + edge_id ∈ main_network_source_edges + allocated = JuMP.value(F[edge_id]) + subnetwork_demands[edge_id][priority_idx] += allocated + end + + user_node_id = edge_id[2] + + if graph[user_node_id].type == :user + allocated = JuMP.value(F[edge_id]) + user_idx = findsorted(user.node_id, user_node_id) + user.allocated[user_idx][priority_idx] = allocated + + # Save allocations to record + push!(record.time, t) + push!(record.allocation_network_id, allocation_model.allocation_network_id) + push!(record.user_node_id, Int(user_node_id)) + push!(record.priority, user.priorities[priority_idx]) + push!(record.demand, user.demand[user_idx][priority_idx](t)) + push!(record.allocated, allocated) + # TODO: This is now the last abstraction before the allocation update, + # should be the average abstraction since the last allocation solve + push!( + record.abstracted, + get_flow(graph, inflow_id(graph, user_node_id), user_node_id, 0), + ) + end + end + + # Write the flows to the subnetworks as allocated flows + # in the allocation object + if is_main_network(allocation_network_id) + for (allocation_network_id, main_network_source_edges) in + zip(allocation_network_ids, main_network_connections) + if is_main_network(allocation_network_id) + continue + end + for edge_id in main_network_source_edges + subnetwork_allocateds[edge_id][priority_idx] = JuMP.value(F[edge_id]) + end + end + end + return nothing +end + +""" +Adjust the source flows. +""" +function adjust_source_capacities!( + allocation_model::AllocationModel, + p::Parameters, + priority_idx::Int; + collect_demands::Bool = false, +)::Nothing + (; problem) = allocation_model + (; graph, allocation) = p + (; allocation_network_id) = allocation_model + (; subnetwork_allocateds) = allocation + edge_ids = graph[].edge_ids[allocation_network_id] + source_constraints = problem[:source] + F = problem[:F] + + main_network_source_edges = get_main_network_connections(p, allocation_network_id) + + 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 priority_idx == 1 + # If the optimization was just started, i.e. sources have to be reset + if edge_id in main_network_source_edges + if collect_demands + # Set the source capacity to effectively unlimited if subnetwork demands are being collected + source_capacity = Inf + else + # Set the source capacity to the value allocated to the subnetwork over this edge + source_capacity = subnetwork_allocateds[edge_id][priority_idx] + end + else + # Reset the source to the current flow from the physical layer. + source_capacity = get_flow(graph, edge_id..., 0) + end + JuMP.set_normalized_rhs( + source_constraints[edge_id], + # It is assumed that the allocation procedure does not have to be differentiated. + source_capacity, + ) + else + # Subtract the allocated flow from the source. + JuMP.set_normalized_rhs( + source_constraints[edge_id], + JuMP.normalized_rhs(source_constraints[edge_id]) - + JuMP.value(F[edge_id]), + ) + end + end + end + return nothing +end + +""" +Set the values of the edge capacities. 2 cases: +- Before the first allocation solve, set the edge capacities to their full capacity; +- Before an allocation solve, subtract the flow used by allocation for the previous priority + from the edge capacities. +""" +function adjust_edge_capacities!( + allocation_model::AllocationModel, + p::Parameters, + priority_idx::Int, +)::Nothing + (; graph) = p + (; problem, capacity, allocation_network_id) = allocation_model + edge_ids = graph[].edge_ids[allocation_network_id] + constraints_capacity = problem[:capacity] + F = problem[:F] + + main_network_source_edges = get_main_network_connections(p, allocation_network_id) + + for edge_id in edge_ids + c = capacity[edge_id...] + + # These edges have no capacity constraints: + # - With infinite capacity + # - Being a source from the main network to a subnetwork + if isinf(c) || edge_id ∈ main_network_source_edges + continue + end + + if priority_idx == 1 + # Before the first allocation solve, set the edge capacities to their full capacity + JuMP.set_normalized_rhs(constraints_capacity[edge_id], c) + else + # Before an allocation solve, subtract the flow used by allocation for the previous priority + # from the edge capacities + JuMP.set_normalized_rhs( + constraints_capacity[edge_id], + JuMP.normalized_rhs(constraints_capacity[edge_id]) - JuMP.value(F[edge_id]), + ) + end + end +end + +function get_basin_capacity( + allocation_model::AllocationModel, + u::ComponentVector, + p::Parameters, + t::Float64, + node_id::NodeID, +)::Float64 + (; graph, basin, allocation_level_control) = p + (; Δt_allocation) = allocation_model + influx = get_flow(graph, node_id, 0.0) + _, basin_idx = id_index(basin.node_id, node_id) + storage_basin = u.storage[basin_idx] + allocation_level_control_node_id = + first(inneighbor_labels_type(graph, node_id, EdgeType.control)) + allocation_level_control_idx = + findsorted(allocation_level_control.node_id, allocation_level_control_node_id) + level_max = allocation_level_control.max_level[allocation_level_control_idx](t) + storage_max = get_storage_from_level(basin, basin_idx, level_max) + + return max(0.0, (storage_basin - storage_max) / Δt_allocation + influx) +end + +function adjust_basin_capacities!( + allocation_model::AllocationModel, + u::ComponentVector, + p::Parameters, + t::Float64, + priority_idx::Int, +)::Nothing + (; problem) = allocation_model + constraints_basin_flow = problem[:basin_flow] + F_basin = problem[:F_basin] + + for node_id in only(constraints_basin_flow.axes) + constraint = constraints_basin_flow[node_id] + if priority_idx == 1 + JuMP.set_normalized_rhs( + constraint, + -get_basin_capacity(allocation_model, u, p, t, node_id), + ) + else + JuMP.set_normalized_rhs( + constraint, + JuMP.normalized_rhs(constraint) - JuMP.value(F_basin[node_id]), + ) + end + end + + return nothing +end + +""" +Save the allocation flows per physical edge. +""" +function save_allocation_flows!( + p::Parameters, + t::Float64, + allocation_model::AllocationModel, + priority::Int, + collect_demands::Bool, +)::Nothing + (; problem, allocation_network_id) = allocation_model + (; allocation, graph) = p + (; record) = allocation + F = problem[:F] + + for allocation_edge in first(F.axes) + flow = JuMP.value(F[allocation_edge]) + edge_metadata = graph[allocation_edge...] + (; node_ids) = edge_metadata + + for i in eachindex(node_ids)[1:(end - 1)] + push!(record.time, t) + push!(record.edge_id, edge_metadata.id) + push!(record.from_node_id, node_ids[i]) + push!(record.to_node_id, node_ids[i + 1]) + push!(record.allocation_network_id, allocation_network_id) + push!(record.priority, priority) + push!(record.flow, flow) + push!(record.collect_demands, collect_demands) + end + end + return nothing +end + +""" +Update the allocation optimization problem for the given subnetwork with the problem state +and flows, solve the allocation problem and assign the results to the users. +""" +function allocate!( + p::Parameters, + allocation_model::AllocationModel, + t::Float64, + u::ComponentVector; + collect_demands::Bool = false, +)::Nothing + (; user, allocation) = p + (; problem, allocation_network_id) = allocation_model + (; priorities) = user + (; subnetwork_demands) = allocation + + main_network_source_edges = get_main_network_connections(p, allocation_network_id) + + if collect_demands + for main_network_connection in keys(subnetwork_demands) + if main_network_connection in main_network_source_edges + subnetwork_demands[main_network_connection] .= 0.0 + end + end + end + + for priority_idx in eachindex(priorities) + adjust_source_capacities!(allocation_model, p, priority_idx; collect_demands) + + # Subtract the flows used by the allocation of the previous priority from the capacities of the edges + # or set edge capacities if priority_idx = 1 + adjust_edge_capacities!(allocation_model, p, priority_idx) + + adjust_basin_capacities!(allocation_model, u, p, t, priority_idx) + + # Set the objective depending on the demands + # A new objective function is set instead of modifying the coefficients + # of an existing objective function because this is not supported for + # quadratic terms: + # https://jump.dev/JuMP.jl/v1.16/manual/objective/#Modify-an-objective-coefficient + set_objective_priority!(allocation_model, p, t, priority_idx) + + # Solve the allocation problem for this priority + JuMP.optimize!(problem) + @debug JuMP.solution_summary(problem) + if JuMP.termination_status(problem) !== JuMP.OPTIMAL + (; allocation_network_id) = allocation_model + priority = priorities[priority_index] + error( + "Allocation of subnetwork $allocation_network_id, priority $priority coudn't find optimal solution.", + ) + end + + # Assign the allocations to the users for this priority + assign_allocations!(allocation_model, p, t, priority_idx; collect_demands) + + # Save the flows over all edges in the subnetwork + save_allocation_flows!( + p, + t, + allocation_model, + priorities[priority_idx], + collect_demands, + ) + end +end From e5b7c88c1f2b2e4a20a76aa1b7e15161109a5dd6 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 7 Feb 2024 17:28:48 +0100 Subject: [PATCH 12/37] Refactor priority handling --- core/src/allocation_init.jl | 4 ++-- core/src/create.jl | 5 +++-- core/src/solve.jl | 13 ++----------- core/src/utils.jl | 15 +++++++++++++++ 4 files changed, 22 insertions(+), 15 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 77a923633..506259dfd 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -1,7 +1,7 @@ """Find the edges from the main network to a subnetwork.""" function find_subnetwork_connections!(p::Parameters)::Nothing - (; allocation, graph, user) = p - n_priorities = length(user.demand[1]) + (; allocation, graph, allocation) = p + n_priorities = length(allocation.priorities) (; subnetwork_demands, subnetwork_allocateds) = allocation for node_id in graph[].node_ids[1] for outflow_id in outflow_ids(graph, node_id) diff --git a/core/src/create.jl b/core/src/create.jl index 5c39a83fd..be7bb2960 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -644,8 +644,8 @@ function User(db::DB, config::Config)::User error("Problems encountered when parsing User static and time node IDs.") end - # All provided priorities - priorities = sort(unique(union(static.priority, time.priority))) + # All priorities used in the model + priorities = get_all_priorities(db, config) active = BitVector() min_level = Float64[] @@ -837,6 +837,7 @@ function Parameters(db::DB, config::Config)::Parameters Int[], AllocationModel[], Vector{Tuple{NodeID, NodeID}}[], + get_all_priorities(db, config), Dict{Tuple{NodeID, NodeID}, Float64}(), Dict{Tuple{NodeID, NodeID}, Float64}(), (; diff --git a/core/src/solve.jl b/core/src/solve.jl index 389185494..ce6e3cfe8 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -35,6 +35,7 @@ struct Allocation allocation_network_ids::Vector{Int} allocation_models::Vector{AllocationModel} main_network_connections::Vector{Vector{Tuple{NodeID, NodeID}}} + priorities::Vector{Int} subnetwork_demands::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} subnetwork_allocateds::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} record::@NamedTuple{ @@ -436,7 +437,6 @@ struct User <: AbstractParameterNode allocated::Vector{Vector{Float64}} return_factor::Vector{Float64} min_level::Vector{Float64} - priorities::Vector{Int} record::@NamedTuple{ time::Vector{Float64}, allocation_network_id::Vector{Int}, @@ -458,16 +458,7 @@ struct User <: AbstractParameterNode record, ) if valid_demand(node_id, demand, priorities) - return new( - node_id, - active, - demand, - allocated, - return_factor, - min_level, - priorities, - record, - ) + return new(node_id, active, demand, allocated, return_factor, min_level, record) else error("Invalid demand") end diff --git a/core/src/utils.jl b/core/src/utils.jl index b90f08237..9212063d5 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -1282,3 +1282,18 @@ end function is_main_network(allocation_network_id::Int)::Bool return allocation_network_id == 1 end + +function get_all_priorities(db::DB, config::Config)::Vector{Int} + priorities = Set{Int}() + + # TODO: Is there a way to automatically grab all tables with a priority column? + for type in [ + UserStaticV1, + UserTimeV1, + AllocationLevelControlStaticV1, + AllocationLevelControlTimeV1, + ] + union!(priorities, load_structvector(db, config, type).priority) + end + return sort(unique(priorities)) +end From 3f9c271053093c66617d349acf9fb88439a59b3c Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 7 Feb 2024 17:32:08 +0100 Subject: [PATCH 13/37] Towards adding terms in objective function for allocation to basins --- core/src/allocation_optim.jl | 58 +++++++++++++++---- .../ribasim_testmodels/allocation.py | 2 +- 2 files changed, 47 insertions(+), 13 deletions(-) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index ca218b76e..24305eeff 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -63,10 +63,11 @@ function set_objective_priority!( priority_idx::Int, )::Nothing (; objective_type, problem, allocation_network_id) = allocation_model - (; graph, user, allocation) = p + (; graph, user, allocation, basin) = p (; demand, node_id) = user (; main_network_connections, subnetwork_demands) = allocation edge_ids = graph[].edge_ids[allocation_network_id] + node_ids = graph[].node_ids[allocation_network_id] F = problem[:F] if objective_type in [:quadratic_absolute, :quadratic_relative] @@ -117,6 +118,14 @@ function set_objective_priority!( end end + # Add + F_basin = problem[:F_basin] + for node_id in only(F_basin.axes) + if priority_idx == findsorted + d = get_basin_demand(allocation_model, u, p, t, node_id) + end + end + new_objective = JuMP.@expression(problem, ex) JuMP.@objective(problem, Min, new_objective) return nothing @@ -165,7 +174,7 @@ function assign_allocations!( push!(record.time, t) push!(record.allocation_network_id, allocation_model.allocation_network_id) push!(record.user_node_id, Int(user_node_id)) - push!(record.priority, user.priorities[priority_idx]) + push!(record.priority, allocation.priorities[priority_idx]) push!(record.demand, user.demand[user_idx][priority_idx](t)) push!(record.allocated, allocated) # TODO: This is now the last abstraction before the allocation update, @@ -290,13 +299,12 @@ function adjust_edge_capacities!( end end -function get_basin_capacity( +function get_basin_data( allocation_model::AllocationModel, - u::ComponentVector, p::Parameters, - t::Float64, + u::ComponentVector, node_id::NodeID, -)::Float64 +) (; graph, basin, allocation_level_control) = p (; Δt_allocation) = allocation_model influx = get_flow(graph, node_id, 0.0) @@ -306,12 +314,39 @@ function get_basin_capacity( first(inneighbor_labels_type(graph, node_id, EdgeType.control)) allocation_level_control_idx = findsorted(allocation_level_control.node_id, allocation_level_control_node_id) - level_max = allocation_level_control.max_level[allocation_level_control_idx](t) - storage_max = get_storage_from_level(basin, basin_idx, level_max) + return storage_basin, Δt_allocation, influx, allocation_level_control_idx, basin_idx +end +function get_basin_capacity( + allocation_model::AllocationModel, + u::ComponentVector, + p::Parameters, + t::Float64, + node_id::NodeID, +)::Float64 + (; allocation_level_control) = p + storage_basin, Δt_allocation, influx, allocation_level_control_idx, basin_idx = + get_basin_data(allocation_model, p, u, node_id) + level_max = allocation_level_control.max_level[allocation_level_control_idx](t) + storage_max = get_storage_from_level(p.basin, basin_idx, level_max) return max(0.0, (storage_basin - storage_max) / Δt_allocation + influx) end +function get_basin_demand( + allocation_model::AllocationModel, + u::ComponentVector, + p::Parameters, + t::Float64, + node_id::NodeID, +)::Float64 + (; allocation_level_control) = p + storage_basin, Δt_allocation, influx, allocation_level_control_idx, basin_idx = + get_basin_data(allocation_model, p, u, node_id) + level_min = allocation_level_control.min_level[allocation_level_control_idx](t) + storage_min = get_storage_from_level(p.basin, basin_idx, level_min) + return max(0.0, (storage_basin - storage_min) / Δt_allocation - influx) +end + function adjust_basin_capacities!( allocation_model::AllocationModel, u::ComponentVector, @@ -386,10 +421,9 @@ function allocate!( u::ComponentVector; collect_demands::Bool = false, )::Nothing - (; user, allocation) = p + (; allocation) = p (; problem, allocation_network_id) = allocation_model - (; priorities) = user - (; subnetwork_demands) = allocation + (; priorities, subnetwork_demands) = allocation main_network_source_edges = get_main_network_connections(p, allocation_network_id) @@ -422,7 +456,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/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 1f3180bc0..9c9ccfcc3 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -1822,7 +1822,7 @@ def allocation_level_control_model(): data={ "node_id": [3], "priority": [2], - "demand": [1.0], + "demand": [1e-3], "return_factor": [0.2], "min_level": [0.2], } From 0edc31b204d741a03b4f5e6477f6a59e4886f821 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 8 Feb 2024 14:19:19 +0100 Subject: [PATCH 14/37] Split basin flow into inflow and outflow --- core/src/allocation_init.jl | 42 +++++++++++++++++++++++++++--------- core/src/allocation_optim.jl | 2 +- 2 files changed, 33 insertions(+), 11 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 506259dfd..9a3f7f98e 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -415,7 +415,10 @@ function add_variables_basin!( node_id for node_id in graph[].node_ids[allocation_network_id] if graph[node_id].type == :basin ] - problem[:F_basin] = JuMP.@variable(problem, F_basin[node_id = node_ids_basin,]) + problem[:F_basin_in] = + JuMP.@variable(problem, F_basin_in[node_id = node_ids_basin,] >= 0.0) + problem[:F_basin_out] = + JuMP.@variable(problem, F_basin_out[node_id = node_ids_basin,] >= 0.0) return nothing end @@ -541,13 +544,25 @@ function outflow_ids_allocation(graph::MetaGraph, node_id::NodeID) return outflow_ids end -function get_basin_flow( +function get_basin_inflow( problem::JuMP.Model, node_id::NodeID, )::Union{JuMP.VariableRef, Float64} - F_basin = problem[:F_basin] - return if node_id in only(F_basin.axes) - F_basin[node_id] + F_basin_in = problem[:F_basin_in] + return if node_id in only(F_basin_in.axes) + F_basin_in[node_id] + else + 0.0 + end +end + +function get_basin_outflow( + problem::JuMP.Model, + node_id::NodeID, +)::Union{JuMP.VariableRef, Float64} + F_basin_out = problem[:F_basin_out] + return if node_id in only(F_basin_out.axes) + F_basin_out[node_id] else 0.0 end @@ -582,11 +597,11 @@ function add_constraints_flow_conservation!( sum([ F[(node_id, outneighbor_id)] for outneighbor_id in outflow_ids_allocation(graph, node_id) - ]) <= + ]) + get_basin_outflow(problem, node_id) <= sum([ F[(inneighbor_id, node_id)] for inneighbor_id in inflow_ids_allocation(graph, node_id) - ]) + get_basin_flow(problem, node_id), + ]) + get_basin_inflow(problem, node_id), base_name = "flow_conservation", ) return nothing @@ -759,12 +774,19 @@ Constraint: flow basin >= - basin capacity """ function add_constraints_basin_flow!(problem::JuMP.Model)::Nothing - F_basin = problem[:F_basin] + F_basin_in = problem[:F_basin_in] + F_basin_out = problem[:F_basin_out] + problem[:basin_flow] = JuMP.@constraint( + problem, + [node_id = only(F_basin.axes)], + F_basin_in[node_id] <= 0.0, + base_name = "basin_inflow" + ) problem[:basin_flow] = JuMP.@constraint( problem, [node_id = only(F_basin.axes)], - F_basin[node_id] >= 0.0, - base_name = "basin_flow" + F_basin_out[node_id] <= 0.0, + base_name = "basin_outflow" ) return nothing end diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 24305eeff..ab1291166 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -363,7 +363,7 @@ function adjust_basin_capacities!( if priority_idx == 1 JuMP.set_normalized_rhs( constraint, - -get_basin_capacity(allocation_model, u, p, t, node_id), + get_basin_capacity(allocation_model, u, p, t, node_id), ) else JuMP.set_normalized_rhs( From e5be0f01006239bd0d8d570923bf7be5727d936d Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 8 Feb 2024 14:30:22 +0100 Subject: [PATCH 15/37] Towards saving basin allocation flows --- core/src/allocation_optim.jl | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index ab1291166..f127b515b 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -377,7 +377,7 @@ function adjust_basin_capacities!( end """ -Save the allocation flows per physical edge. +Save the allocation flows per basin physical edge. """ function save_allocation_flows!( p::Parameters, @@ -391,6 +391,7 @@ function save_allocation_flows!( (; record) = allocation F = problem[:F] + # Edge flows for allocation_edge in first(F.axes) flow = JuMP.value(F[allocation_edge]) edge_metadata = graph[allocation_edge...] @@ -407,6 +408,22 @@ function save_allocation_flows!( push!(record.collect_demands, collect_demands) end end + + # Basin flows + for node_id in graph[].node_id[allocation_network_id] + if graph[node_id].type == :basin + flow = ... + push!(record.time, t) + push!(record.edge_id, 0) + push!(record.from_node_id, node_id) + push!(record.to_node_id, node_id) + push!(record.allocation_network_id, allocation_network_id) + push!(record.priority, priority) + push!(record.flow, flow) + push!(record.collect_demands, collect_demands) + end + end + return nothing end From bfcc07d66fe3acb6fcf721dae78818af91a1f245 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 8 Feb 2024 16:54:58 +0100 Subject: [PATCH 16/37] Towards adding flows into basins to objective function --- core/src/allocation_init.jl | 11 +-- core/src/allocation_optim.jl | 146 ++++++++++++++++++++++++++--------- core/src/utils.jl | 11 +++ 3 files changed, 121 insertions(+), 47 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 9a3f7f98e..66a352f5a 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -582,7 +582,6 @@ function add_constraints_flow_conservation!( )::Nothing (; graph) = p F = problem[:F] - F_basin = problem[:F_basin] node_ids = graph[].node_ids[allocation_network_id] node_ids_conservation = [node_id for node_id in node_ids if graph[node_id].type == :basin] @@ -776,15 +775,9 @@ flow basin >= - basin capacity function add_constraints_basin_flow!(problem::JuMP.Model)::Nothing F_basin_in = problem[:F_basin_in] F_basin_out = problem[:F_basin_out] - problem[:basin_flow] = JuMP.@constraint( + problem[:basin_outflow] = JuMP.@constraint( problem, - [node_id = only(F_basin.axes)], - F_basin_in[node_id] <= 0.0, - base_name = "basin_inflow" - ) - problem[:basin_flow] = JuMP.@constraint( - problem, - [node_id = only(F_basin.axes)], + [node_id = only(F_basin_out.axes)], F_basin_out[node_id] <= 0.0, base_name = "basin_outflow" ) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index f127b515b..bbc09b1cf 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -1,23 +1,16 @@ -""" -Add a term to the expression of the objective function corresponding to -the demand of a user. -""" -function add_user_term!( +function add_objective_term!( ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, - edge::Tuple{NodeID, NodeID}, - objective_type::Symbol, + problem::JuMP.Model, + F_variable::JuMP.VariableRef, demand::Float64, - model::AllocationModel, + objective_type::Symbol; + constraint_abs_positive::Union{JuMP.ConstraintRef, Nothing} = nothing, + constraint_abs_negative::Union{JuMP.ConstraintRef, Nothing} = nothing, )::Nothing - (; problem) = model - F = problem[:F] - F_edge = F[edge] - node_id_user = edge[2] - if objective_type == :quadratic_absolute # Objective function ∑ (F - d)^2 - JuMP.add_to_expression!(ex, 1, F_edge, F_edge) - JuMP.add_to_expression!(ex, -2 * demand, F_edge) + JuMP.add_to_expression!(ex, 1, F_variable, F_variable) + JuMP.add_to_expression!(ex, -2 * demand, F_variable) JuMP.add_to_expression!(ex, demand^2) elseif objective_type == :quadratic_relative @@ -25,24 +18,24 @@ function add_user_term!( if demand ≈ 0 return nothing end - JuMP.add_to_expression!(ex, 1.0 / demand^2, F_edge, F_edge) - JuMP.add_to_expression!(ex, -2.0 / demand, F_edge) + JuMP.add_to_expression!(ex, 1.0 / demand^2, F_variable, F_variable) + JuMP.add_to_expression!(ex, -2.0 / demand, F_variable) JuMP.add_to_expression!(ex, 1.0) elseif objective_type == :linear_absolute # Objective function ∑ |F - d| - JuMP.set_normalized_rhs(problem[:abs_positive][node_id_user], -demand) - JuMP.set_normalized_rhs(problem[:abs_negative][node_id_user], demand) + JuMP.set_normalized_rhs(constraint_abs_positive, -demand) + JuMP.set_normalized_rhs(constraint_abs_negative, demand) elseif objective_type == :linear_relative # Objective function ∑ |1 - F/d| JuMP.set_normalized_coefficient( - problem[:abs_positive][node_id_user], + constraint_abs_positive, F_edge, iszero(demand) ? 0 : 1 / demand, ) JuMP.set_normalized_coefficient( - problem[:abs_negative][node_id_user], + constraint_abs_negative, F_edge, iszero(demand) ? 0 : -1 / demand, ) @@ -52,6 +45,74 @@ function add_user_term!( return nothing end +""" +Add a term to the expression of the objective function corresponding to +the demand of a user. +""" +function add_user_term!( + ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, + edge::Tuple{NodeID, NodeID}, + objective_type::Symbol, + demand::Float64, + problem::JuMP.Model, +)::Nothing + F = problem[:F] + F_edge = F[edge] + node_id_user = edge[2] + + if objective_type in [:linear_absolute, :linear_relative] + constraint_abs_positive = problem[:abs_positive][node_id_user] + constraint_abs_negative = problem[:abs_negative][node_id_user] + else + constraint_abs_positive = nothing + constraint_abs_negative = nothing + end + + add_objective_term!( + ex, + problem, + F_edge, + demand, + objective_type; + 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!( + ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, + problem::JuMP.Model, + demand::Float64, + objective_type::Symbol, + node_id::NodeID, +)::Nothing + F_basin_in = problem[:F_basin_in] + F_basin = F_basin_in[node_id] + + if objective_type in [:linear_absolute, :linear_relative] + constraint_abs_positive = problem[:abs_positive][node_id_user] + constraint_abs_negative = problem[:abs_negative][node_id_user] + else + constraint_abs_positive = nothing + constraint_abs_negative = nothing + end + + add_objective_term!( + ex, + problem, + F_basin, + demand, + objective_type; + constraint_abs_positive, + constraint_abs_negative, + ) + return nothing +end + """ Set the objective for the given priority. For an objective with absolute values this also involves adjusting constraints. @@ -59,13 +120,14 @@ For an objective with absolute values this also involves adjusting constraints. function set_objective_priority!( allocation_model::AllocationModel, p::Parameters, + u::ComponentVector, t::Float64, priority_idx::Int, )::Nothing (; objective_type, problem, allocation_network_id) = allocation_model (; graph, user, allocation, basin) = p (; demand, node_id) = user - (; main_network_connections, subnetwork_demands) = allocation + (; main_network_connections, subnetwork_demands, priorities) = allocation edge_ids = graph[].edge_ids[allocation_network_id] node_ids = graph[].node_ids[allocation_network_id] @@ -85,7 +147,7 @@ function set_objective_priority!( 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) + add_user_term!(ex, connection, objective_type, d, problem) end end end @@ -100,7 +162,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) - add_user_term!(ex, edge_id, objective_type, d, allocation_model) + add_user_term!(ex, edge_id, objective_type, d, problem) end # Add flow cost @@ -118,12 +180,18 @@ function set_objective_priority!( end end - # Add - F_basin = problem[:F_basin] - for node_id in only(F_basin.axes) - if priority_idx == findsorted - d = get_basin_demand(allocation_model, u, p, t, node_id) - end + # Terms for basins + F_basin_in = problem[:F_basin_in] + F_basin_out = problem[:F_basin_out] + for node_id in only(F_basin_in.axes) + priority_basin = get_basin_priority(p, node_id) + priority_now = priorities[priority_idx] + + d = + priority_basin == priority_now ? + get_basin_demand(allocation_model, u, p, t, node_id) : 0.0 + + #add_objective_term() end new_objective = JuMP.@expression(problem, ex) @@ -355,11 +423,11 @@ function adjust_basin_capacities!( priority_idx::Int, )::Nothing (; problem) = allocation_model - constraints_basin_flow = problem[:basin_flow] - F_basin = problem[:F_basin] + constraints_outflow = problem[:basin_outflow] + F_basin_out = problem[:F_basin_out] - for node_id in only(constraints_basin_flow.axes) - constraint = constraints_basin_flow[node_id] + for node_id in only(constraints_outflow.axes) + constraint = constraints_outflow[node_id] if priority_idx == 1 JuMP.set_normalized_rhs( constraint, @@ -368,7 +436,7 @@ function adjust_basin_capacities!( else JuMP.set_normalized_rhs( constraint, - JuMP.normalized_rhs(constraint) - JuMP.value(F_basin[node_id]), + JuMP.normalized_rhs(constraint) - JuMP.value(F_basin_out[node_id]), ) end end @@ -390,6 +458,8 @@ function save_allocation_flows!( (; allocation, graph) = p (; record) = allocation F = problem[:F] + F_basin_in = problem[:F_basin_in] + F_basin_out = problem[:F_basin_out] # Edge flows for allocation_edge in first(F.axes) @@ -410,9 +480,9 @@ function save_allocation_flows!( end # Basin flows - for node_id in graph[].node_id[allocation_network_id] + for node_id in graph[].node_ids[allocation_network_id] if graph[node_id].type == :basin - flow = ... + flow = JuMP.value(F_basin_out[node_id]) - JuMP.value(F_basin_in[node_id]) push!(record.time, t) push!(record.edge_id, 0) push!(record.from_node_id, node_id) @@ -466,7 +536,7 @@ function allocate!( # of an existing objective function because this is not supported for # quadratic terms: # https://jump.dev/JuMP.jl/v1.16/manual/objective/#Modify-an-objective-coefficient - set_objective_priority!(allocation_model, p, t, priority_idx) + set_objective_priority!(allocation_model, p, u, t, priority_idx) # Solve the allocation problem for this priority JuMP.optimize!(problem) diff --git a/core/src/utils.jl b/core/src/utils.jl index 9212063d5..84fadff59 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -1297,3 +1297,14 @@ function get_all_priorities(db::DB, config::Config)::Vector{Int} end return sort(unique(priorities)) end + +function get_basin_priority(p::Parameters, node_id::NodeID)::Int + (; graph, allocation_level_control) = p + inneighbors_control = inneighbor_labels_type(graph, node_id, EdgeType.control) + if isempty(inneighbors_control) + return 0 + else + idx = findsorted(allocation_level_control.node_id, first(inneighbors_control)) + return allocation_level_control.priority[idx] + end +end From a560d43861142eedf3af2c457a9a958ee97f6557 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 9 Feb 2024 09:48:51 +0100 Subject: [PATCH 17/37] Add basin absolute value variables --- core/src/allocation_init.jl | 34 +++++++++++++++++++++++----------- core/src/allocation_optim.jl | 34 +++++++++++++++++++++++----------- core/test/allocation_test.jl | 30 ++++++++++++++++-------------- 3 files changed, 62 insertions(+), 36 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 66a352f5a..97c12c7c1 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -438,7 +438,17 @@ function add_variables_absolute_value!( (; main_network_connections) = allocation if startswith(config.allocation.objective_type, "linear") node_ids = graph[].node_ids[allocation_network_id] - node_ids_user = [node_id for node_id in node_ids if graph[node_id].type == :user] + node_ids_user = NodeID[] + node_ids_basin = NodeID[] + + for node_id in node_ids + type = graph[node_id].type + if type == :user + push!(node_ids_user, node_id) + elseif type == :basin + push!(node_ids_basin, node_id) + end + end # For the main network, connections to subnetworks are treated as users if is_main_network(allocation_network_id) @@ -449,7 +459,9 @@ function add_variables_absolute_value!( end end - problem[:F_abs] = JuMP.@variable(problem, F_abs[node_id = node_ids_user]) + problem[:F_abs_user] = JuMP.@variable(problem, F_abs_user[node_id = node_ids_user]) + problem[:F_abs_basin] = + JuMP.@variable(problem, F_abs_basin[node_id = node_ids_basin]) end return nothing end @@ -671,40 +683,40 @@ function add_constraints_absolute_value!( node_id_user in node_ids_user ) F = problem[:F] - F_abs = problem[:F_abs] + F_abs_user = problem[:F_abs_user] d = 2.0 if config.allocation.objective_type == "linear_absolute" - # 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) + # These constraints together make sure that F_abs_user acts as the absolute + # value F_abs_user = |x| where x = F-d (here for example d = 2) problem[:abs_positive] = JuMP.@constraint( problem, [node_id_user = node_ids_user], - F_abs[node_id_user] >= + F_abs_user[node_id_user] >= (F[(node_ids_user_inflow[node_id_user], node_id_user)] - d), base_name = "abs_positive" ) problem[:abs_negative] = JuMP.@constraint( problem, [node_id_user = node_ids_user], - F_abs[node_id_user] >= + F_abs_user[node_id_user] >= -(F[(node_ids_user_inflow[node_id_user], node_id_user)] - d), base_name = "abs_negative" ) elseif config.allocation.objective_type == "linear_relative" - # These constraints together make sure that F_abs acts as the absolute - # value F_abs = |x| where x = 1-F/d (here for example d = 2) + # These constraints together make sure that F_abs_user acts as the absolute + # value F_abs_user = |x| where x = 1-F/d (here for example d = 2) problem[:abs_positive] = JuMP.@constraint( problem, [node_id_user = node_ids_user], - F_abs[node_id_user] >= + F_abs_user[node_id_user] >= (1 - F[(node_ids_user_inflow[node_id_user], node_id_user)] / d), base_name = "abs_positive" ) problem[:abs_negative] = JuMP.@constraint( problem, [node_id_user = node_ids_user], - F_abs[node_id_user] >= + F_abs_user[node_id_user] >= -(1 - F[(node_ids_user_inflow[node_id_user], node_id_user)] / d), base_name = "abs_negative" ) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index bbc09b1cf..968962dac 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -136,7 +136,7 @@ function set_objective_priority!( ex = JuMP.QuadExpr() elseif objective_type in [:linear_absolute, :linear_relative] ex = JuMP.AffExpr() - ex += sum(problem[:F_abs]) + ex += sum(problem[:F_abs_user]) end demand_max = 0.0 @@ -378,10 +378,14 @@ function get_basin_data( influx = get_flow(graph, node_id, 0.0) _, basin_idx = id_index(basin.node_id, node_id) storage_basin = u.storage[basin_idx] - allocation_level_control_node_id = - first(inneighbor_labels_type(graph, node_id, EdgeType.control)) - allocation_level_control_idx = - findsorted(allocation_level_control.node_id, allocation_level_control_node_id) + control_inneighbors = inneighbor_labels_type(graph, node_id, EdgeType.control) + if isempty(control_inneighbors) + allocation_level_control_idx = 0 + else + allocation_level_control_node_id = first(control_inneighbors) + allocation_level_control_idx = + findsorted(allocation_level_control.node_id, allocation_level_control_node_id) + end return storage_basin, Δt_allocation, influx, allocation_level_control_idx, basin_idx end @@ -395,9 +399,13 @@ function get_basin_capacity( (; allocation_level_control) = p storage_basin, Δt_allocation, influx, allocation_level_control_idx, basin_idx = get_basin_data(allocation_model, p, u, node_id) - level_max = allocation_level_control.max_level[allocation_level_control_idx](t) - storage_max = get_storage_from_level(p.basin, basin_idx, level_max) - return max(0.0, (storage_basin - storage_max) / Δt_allocation + influx) + if iszero(allocation_level_control_idx) + return 0.0 + else + level_max = allocation_level_control.max_level[allocation_level_control_idx](t) + storage_max = get_storage_from_level(p.basin, basin_idx, level_max) + return max(0.0, (storage_basin - storage_max) / Δt_allocation + influx) + end end function get_basin_demand( @@ -410,9 +418,13 @@ function get_basin_demand( (; allocation_level_control) = p storage_basin, Δt_allocation, influx, allocation_level_control_idx, basin_idx = get_basin_data(allocation_model, p, u, node_id) - level_min = allocation_level_control.min_level[allocation_level_control_idx](t) - storage_min = get_storage_from_level(p.basin, basin_idx, level_min) - return max(0.0, (storage_basin - storage_min) / Δt_allocation - influx) + if iszero(allocation_level_control_idx) + return 0.0 + else + level_min = allocation_level_control.min_level[allocation_level_control_idx](t) + storage_min = get_storage_from_level(p.basin, basin_idx, level_min) + return max(0.0, (storage_basin - storage_min) / Δt_allocation - influx) + end end function adjust_basin_capacities!( diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 35632e181..6edd008bd 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -1,5 +1,6 @@ @testitem "Allocation solve" begin using Ribasim: NodeID + using ComponentArrays: ComponentVector import SQLite import JuMP @@ -23,9 +24,10 @@ end end + u = ComponentVector(; storage = zeros(length(p.basin.node_id))) Ribasim.set_flow!(graph, NodeID(1), NodeID(2), 4.5) # Source flow allocation_model = p.allocation.allocation_models[1] - Ribasim.allocate!(p, allocation_model, 0.0) + Ribasim.allocate!(p, allocation_model, 0.0, u) F = allocation_model.problem[:F] @test JuMP.value(F[(NodeID(2), NodeID(6))]) ≈ 0.0 @@ -90,12 +92,12 @@ end problem = model.integrator.p.allocation.allocation_models[1].problem objective = JuMP.objective_function(problem) @test objective isa JuMP.AffExpr # Affine expression - @test :F_abs in keys(problem.obj_dict) + @test :F_abs_user in keys(problem.obj_dict) F = problem[:F] - F_abs = problem[:F_abs] + F_abs_user = problem[:F_abs_user] - @test objective.terms[F_abs[NodeID(5)]] == 1.0 - @test objective.terms[F_abs[NodeID(6)]] == 1.0 + @test objective.terms[F_abs_user[NodeID(5)]] == 1.0 + @test objective.terms[F_abs_user[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 @@ -107,12 +109,12 @@ end problem = model.integrator.p.allocation.allocation_models[1].problem objective = JuMP.objective_function(problem) @test objective isa JuMP.AffExpr # Affine expression - @test :F_abs in keys(problem.obj_dict) + @test :F_abs_user in keys(problem.obj_dict) F = problem[:F] - F_abs = problem[:F_abs] + F_abs_user = problem[:F_abs_user] - @test objective.terms[F_abs[NodeID(5)]] == 1.0 - @test objective.terms[F_abs[NodeID(6)]] == 1.0 + @test objective.terms[F_abs_user[NodeID(5)]] == 1.0 + @test objective.terms[F_abs_user[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 @@ -208,7 +210,7 @@ end # support absolute value expressions in the objective function allocation_model_main_network = Ribasim.get_allocation_model(p, 1) problem = allocation_model_main_network.problem - @test problem[:F_abs].axes[1] == NodeID[11, 24, 38] + @test problem[:F_abs_user].axes[1] == NodeID[11, 24, 38] @test problem[:abs_positive].axes[1] == NodeID[11, 24, 38] @test problem[:abs_negative].axes[1] == NodeID[11, 24, 38] @@ -260,10 +262,10 @@ end # Main network objective function objective = JuMP.objective_function(problem) objective_variables = keys(objective.terms) - F_abs = problem[:F_abs] - @test F_abs[NodeID(11)] ∈ objective_variables - @test F_abs[NodeID(24)] ∈ objective_variables - @test F_abs[NodeID(38)] ∈ objective_variables + F_abs_user = problem[:F_abs_user] + @test F_abs_user[NodeID(11)] ∈ objective_variables + @test F_abs_user[NodeID(24)] ∈ objective_variables + @test F_abs_user[NodeID(38)] ∈ objective_variables # Running full allocation algorithm Ribasim.set_flow!(graph, NodeID(1), NodeID(2), 4.5) From 22d73bbd9aa48c587f0215bba6adddc1433b2eea Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 9 Feb 2024 12:59:59 +0100 Subject: [PATCH 18/37] Add constraints on basin absolute value variables --- core/src/allocation_init.jl | 144 ++++++++++++++++++++++------------- core/src/allocation_optim.jl | 4 +- 2 files changed, 93 insertions(+), 55 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 6f1378d78..ec581538e 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -623,6 +623,56 @@ 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, + objective_type::String, + variable_type::String, +)::Nothing + # Example demand + d = 2.0 + + node_ids = only(F_abs.axes) + + if objective_type == "linear_absolute" + # 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 + ) + elseif objective_type == "linear_relative" + # These constraints together make sure that F_abs_user acts as the absolute + # value F_abs_user = |x| where x = 1-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] >= (1 - 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_user] >= -(1 - flow_per_node[node_id] / d), + base_name = base_name + ) + end + return nothing +end + +function add_constraints_absolute_value_user!( problem::JuMP.Model, p::Parameters, allocation_network_id::Int, @@ -633,61 +683,48 @@ function add_constraints_absolute_value!( objective_type = config.allocation.objective_type if startswith(objective_type, "linear") - node_ids = graph[].node_ids[allocation_network_id] - node_ids_user = [node_id for node_id in node_ids if graph[node_id].type == :user] + F = problem[:F] + F_abs_user = problem[:F_abs_user] - # For the main network, connections to subnetworks are treated as users - if is_main_network(allocation_network_id) - for connections_subnetwork in main_network_connections - for connection in connections_subnetwork - push!(node_ids_user, connection[2]) - end - end - end + flow_per_node = Dict( + node_id => F[(only(inflow_ids_allocation(graph, node_id)), node_id)] for + node_id in only(F_abs_user.axes) + ) - node_ids_user_inflow = Dict( - node_id_user => only(inflow_ids_allocation(graph, node_id_user)) for - node_id_user in node_ids_user + add_constraints_absolute_value!( + problem, + flow_per_node, + F_abs_user, + objective_type, + "user", + ) + end + return nothing +end + +function add_constraints_absolute_value_basin!( + problem::JuMP.Model, + p::Parameters, + allocation_network_id::Int, + config::Config, +)::Nothing + (; graph, allocation) = p + (; main_network_connections) = allocation + + objective_type = config.allocation.objective_type + if startswith(objective_type, "linear") + F_basin_in = problem[:F_basin_in] + F_abs_basin = problem[:F_abs_basin] + flow_per_node = + Dict(node_id => F_basin_in[node_id] for node_id in only(F_abs_basin.axes)) + + add_constraints_absolute_value!( + problem, + flow_per_node, + F_abs_basin, + objective_type, + "basin", ) - F = problem[:F] - F_abs_user = problem[:F_abs_user] - d = 2.0 - - if config.allocation.objective_type == "linear_absolute" - # These constraints together make sure that F_abs_user acts as the absolute - # value F_abs_user = |x| where x = F-d (here for example d = 2) - problem[:abs_positive] = JuMP.@constraint( - problem, - [node_id_user = node_ids_user], - F_abs_user[node_id_user] >= - (F[(node_ids_user_inflow[node_id_user], node_id_user)] - d), - base_name = "abs_positive" - ) - problem[:abs_negative] = JuMP.@constraint( - problem, - [node_id_user = node_ids_user], - F_abs_user[node_id_user] >= - -(F[(node_ids_user_inflow[node_id_user], node_id_user)] - d), - base_name = "abs_negative" - ) - elseif config.allocation.objective_type == "linear_relative" - # These constraints together make sure that F_abs_user acts as the absolute - # value F_abs_user = |x| where x = 1-F/d (here for example d = 2) - problem[:abs_positive] = JuMP.@constraint( - problem, - [node_id_user = node_ids_user], - F_abs_user[node_id_user] >= - (1 - F[(node_ids_user_inflow[node_id_user], node_id_user)] / d), - base_name = "abs_positive" - ) - problem[:abs_negative] = JuMP.@constraint( - problem, - [node_id_user = node_ids_user], - F_abs_user[node_id_user] >= - -(1 - F[(node_ids_user_inflow[node_id_user], node_id_user)] / d), - base_name = "abs_negative" - ) - end end return nothing end @@ -785,7 +822,8 @@ function allocation_problem( add_constraints_source!(problem, p, allocation_network_id) add_constraints_flow_conservation!(problem, p, allocation_network_id) add_constraints_user_returnflow!(problem, p, allocation_network_id) - add_constraints_absolute_value!(problem, p, allocation_network_id, config) + add_constraints_absolute_value_user!(problem, p, allocation_network_id, config) + add_constraints_absolute_value_basin!(problem, p, allocation_network_id, config) add_constraints_fractional_flow!(problem, p, allocation_network_id) add_constraints_basin_flow!(problem) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 1fd72fe2f..40a0425a9 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -61,8 +61,8 @@ function add_user_term!( node_id_user = edge[2] if objective_type in [:linear_absolute, :linear_relative] - constraint_abs_positive = problem[:abs_positive][node_id_user] - constraint_abs_negative = problem[:abs_negative][node_id_user] + constraint_abs_positive = problem[:abs_positive_user][node_id_user] + constraint_abs_negative = problem[:abs_negative_user][node_id_user] else constraint_abs_positive = nothing constraint_abs_negative = nothing From de7cd836afc249ecd5b8a4a134b3344d91b32e5d Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 9 Feb 2024 14:18:47 +0100 Subject: [PATCH 19/37] Finish adding basin flow terms to objective function --- core/src/allocation_optim.jl | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 40a0425a9..687392422 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -1,6 +1,5 @@ function add_objective_term!( ex::Union{JuMP.QuadExpr, JuMP.AffExpr}, - problem::JuMP.Model, F_variable::JuMP.VariableRef, demand::Float64, objective_type::Symbol; @@ -70,7 +69,6 @@ function add_user_term!( add_objective_term!( ex, - problem, F_edge, demand, objective_type; @@ -94,8 +92,8 @@ function add_basin_term!( F_basin = F_basin_in[node_id] if objective_type in [:linear_absolute, :linear_relative] - constraint_abs_positive = problem[:abs_positive][node_id_user] - constraint_abs_negative = problem[:abs_negative][node_id_user] + constraint_abs_positive = problem[:abs_positive_basin][node_id] + constraint_abs_negative = problem[:abs_negative_basin][node_id] else constraint_abs_positive = nothing constraint_abs_negative = nothing @@ -103,7 +101,6 @@ function add_basin_term!( add_objective_term!( ex, - problem, F_basin, demand, objective_type; @@ -125,11 +122,10 @@ function set_objective_priority!( priority_idx::Int, )::Nothing (; objective_type, problem, allocation_network_id) = allocation_model - (; graph, user, allocation, basin) = p - (; demand, demand_itp, demand_from_timeseries, node_id) = user + (; graph, user, allocation) = p + (; demand_itp, demand_from_timeseries, node_id) = user (; main_network_connections, subnetwork_demands, priorities) = allocation edge_ids = graph[].edge_ids[allocation_network_id] - node_ids = graph[].node_ids[allocation_network_id] F = problem[:F] if objective_type in [:quadratic_absolute, :quadratic_relative] @@ -137,6 +133,7 @@ function set_objective_priority!( elseif objective_type in [:linear_absolute, :linear_relative] ex = JuMP.AffExpr() ex += sum(problem[:F_abs_user]) + ex += sum(problem[:F_abs_basin]) end demand_max = 0.0 @@ -187,16 +184,13 @@ function set_objective_priority!( # Terms for basins F_basin_in = problem[:F_basin_in] - F_basin_out = problem[:F_basin_out] for node_id in only(F_basin_in.axes) priority_basin = get_basin_priority(p, node_id) priority_now = priorities[priority_idx] - d = priority_basin == priority_now ? get_basin_demand(allocation_model, u, p, t, node_id) : 0.0 - - #add_objective_term() + add_basin_term!(ex, problem, d, objective_type, node_id) end new_objective = JuMP.@expression(problem, ex) From 59508a0eae1e6572b6663cd774f34c6258b9c789 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 9 Feb 2024 14:36:21 +0100 Subject: [PATCH 20/37] Fix tests --- core/src/allocation_init.jl | 2 +- core/src/allocation_optim.jl | 4 ++-- core/test/allocation_test.jl | 14 ++++++++------ 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index ec581538e..489465ae9 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -665,7 +665,7 @@ function add_constraints_absolute_value!( problem[Symbol(base_name)] = JuMP.@constraint( problem, [node_id = node_ids], - F_abs[node_id_user] >= -(1 - flow_per_node[node_id] / d), + F_abs[node_id] >= -(1 - flow_per_node[node_id] / d), base_name = base_name ) end diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 687392422..301501266 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -30,12 +30,12 @@ function add_objective_term!( # Objective function ∑ |1 - F/d| JuMP.set_normalized_coefficient( constraint_abs_positive, - F_edge, + F_variable, iszero(demand) ? 0 : 1 / demand, ) JuMP.set_normalized_coefficient( constraint_abs_negative, - F_edge, + F_variable, iszero(demand) ? 0 : -1 / demand, ) else diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 3497c9d9d..d29ee1f8f 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -217,8 +217,8 @@ end allocation_model_main_network = Ribasim.get_allocation_model(p, 1) problem = allocation_model_main_network.problem @test problem[:F_abs_user].axes[1] == NodeID[11, 24, 38] - @test problem[:abs_positive].axes[1] == NodeID[11, 24, 38] - @test problem[:abs_negative].axes[1] == NodeID[11, 24, 38] + @test problem[:abs_positive_user].axes[1] == NodeID[11, 24, 38] + @test problem[:abs_negative_user].axes[1] == NodeID[11, 24, 38] # In each subnetwork, the connection from the main network to the subnetwork is # interpreted as a source @@ -233,6 +233,7 @@ end @testitem "allocation with main network optimization problem" begin using SQLite using Ribasim: NodeID + using ComponentArrays: ComponentVector using JuMP toml_path = normpath( @@ -246,13 +247,14 @@ end p = Ribasim.Parameters(db, cfg) close(db) - (; allocation, user, graph) = p + (; allocation, user, graph, basin) = p (; allocation_models, subnetwork_demands, subnetwork_allocateds) = allocation t = 0.0 # Collecting demands + u = ComponentVector(; storage = zeros(length(basin.node_id))) for allocation_model in allocation_models[2:end] - Ribasim.allocate!(p, allocation_model, t; collect_demands = true) + Ribasim.allocate!(p, allocation_model, t, u; collect_demands = true) end @test subnetwork_demands[(NodeID(2), NodeID(11))] ≈ [4.0, 4.0, 0.0] @@ -263,7 +265,7 @@ end # containing subnetworks as users allocation_model = allocation_models[1] (; problem) = allocation_model - Ribasim.allocate!(p, allocation_model, t) + Ribasim.allocate!(p, allocation_model, t, u) # Main network objective function objective = JuMP.objective_function(problem) @@ -275,7 +277,7 @@ end # Running full allocation algorithm Ribasim.set_flow!(graph, NodeID(1), NodeID(2), 4.5) - Ribasim.update_allocation!((; p, t)) + Ribasim.update_allocation!((; p, t, u)) @test subnetwork_allocateds[NodeID(2), NodeID(11)] ≈ [4.0, 0.49766666, 0.0] @test subnetwork_allocateds[NodeID(6), NodeID(24)] ≈ [0.00133333333, 0.0, 0.0] From 1e27b36c073ff3c4b3bd02683ae8e9998049e8e8 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 9 Feb 2024 14:49:39 +0100 Subject: [PATCH 21/37] Fix documentation optimization problem example --- docs/core/allocation.qmd | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index afed1f0f5..98c1916a8 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -234,9 +234,11 @@ The following is an example of an optimization problem for the example shown [he # | code-fold: true using Ribasim using SQLite +using ComponentArrays: ComponentVector toml_path = normpath(@__DIR__, "../../generated_testmodels/allocation_example/ribasim.toml") p = Ribasim.Model(toml_path).integrator.p +u = ComponentVector(; storage = zeros(length(p.basin.node_id))) allocation_model = p.allocation.allocation_models[1] t = 0.0 @@ -246,7 +248,7 @@ Ribasim.set_flow!(p.graph, Ribasim.NodeID(1), Ribasim.NodeID(2), 1.0) Ribasim.adjust_source_capacities!(allocation_model, p, priority_idx) Ribasim.adjust_edge_capacities!(allocation_model, p, priority_idx) -Ribasim.set_objective_priority!(allocation_model, p, t, priority_idx) +Ribasim.set_objective_priority!(allocation_model, p, u, t, priority_idx) println(p.allocation.allocation_models[1].problem) ``` From da374b0489370b1df36f5a9f8660dcb55895b894 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 9 Feb 2024 15:16:39 +0100 Subject: [PATCH 22/37] Add docstrings --- core/src/allocation_init.jl | 26 ++++++++++++-------------- core/src/allocation_optim.jl | 32 ++++++++++++++++++++++++++++++++ core/src/parameter.jl | 6 ++++++ 3 files changed, 50 insertions(+), 14 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 489465ae9..e7c6c47df 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -672,14 +672,16 @@ function add_constraints_absolute_value!( return nothing end +""" +Add constraints so that variables F_abs_user act as the +absolute value of the expression comparing flow to an user to its demand. +""" function add_constraints_absolute_value_user!( problem::JuMP.Model, p::Parameters, - allocation_network_id::Int, config::Config, )::Nothing - (; graph, allocation) = p - (; main_network_connections) = allocation + (; graph) = p objective_type = config.allocation.objective_type if startswith(objective_type, "linear") @@ -702,15 +704,11 @@ function add_constraints_absolute_value_user!( return nothing end -function add_constraints_absolute_value_basin!( - problem::JuMP.Model, - p::Parameters, - allocation_network_id::Int, - config::Config, -)::Nothing - (; graph, allocation) = p - (; main_network_connections) = allocation - +""" +Add constraints so that variables F_abs_basin act as the +absolute value of the expression comparing flow to a basin to its demand. +""" +function add_constraints_absolute_value_basin!(problem::JuMP.Model, config::Config)::Nothing objective_type = config.allocation.objective_type if startswith(objective_type, "linear") F_basin_in = problem[:F_basin_in] @@ -822,8 +820,8 @@ function allocation_problem( add_constraints_source!(problem, p, allocation_network_id) add_constraints_flow_conservation!(problem, p, allocation_network_id) add_constraints_user_returnflow!(problem, p, allocation_network_id) - add_constraints_absolute_value_user!(problem, p, allocation_network_id, config) - add_constraints_absolute_value_basin!(problem, p, allocation_network_id, config) + add_constraints_absolute_value_user!(problem, p, config) + add_constraints_absolute_value_basin!(problem, config) add_constraints_fractional_flow!(problem, p, allocation_network_id) add_constraints_basin_flow!(problem) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 301501266..37d4d4426 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -1,3 +1,7 @@ +""" +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::Union{JuMP.QuadExpr, JuMP.AffExpr}, F_variable::JuMP.VariableRef, @@ -366,6 +370,15 @@ function adjust_edge_capacities!( end end +""" +Get several variables associated with a basin: +- Its current storage +- The allocation update interval +- The influx (sum of instantaneous vertical fluxes of the basin) +- The index of the connected allocation_level_control node (0 if such a + node does not exist) +- The index of the basin +""" function get_basin_data( allocation_model::AllocationModel, p::Parameters, @@ -388,6 +401,13 @@ function get_basin_data( return storage_basin, Δt_allocation, influx, allocation_level_control_idx, basin_idx end +""" +Get the capacity of the basin, i.e. the maximum +flow that can be abstracted from the basin if it is in a +state of surplus storage (0 if no reference levels are provided by +a allocation_level_control node). +Storages are converted to flows by dividing by the allocation timestep. +""" function get_basin_capacity( allocation_model::AllocationModel, u::ComponentVector, @@ -407,6 +427,12 @@ function get_basin_capacity( end end +""" +Get the demand of the basin, i.e. how large a flow the +basin needs to get to its minimum target level (0 if no +reference levels are provided by a allocation_level_control node). +Storages are converted to flows by dividing by the allocation timestep. +""" function get_basin_demand( allocation_model::AllocationModel, u::ComponentVector, @@ -426,6 +452,12 @@ function get_basin_demand( end end +""" +Set the values of the basin outflows. 2 cases: +- Before the first allocation solve, set the capacities to their full capacity if there is surplus storage; +- Before an allocation solve, subtract the flow used by allocation for the previous priority + from the capacities. +""" function adjust_basin_capacities!( allocation_model::AllocationModel, u::ComponentVector, diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 0861054fc..d9ef5bd3e 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -496,6 +496,12 @@ struct User <: AbstractParameterNode end end +""" +node_id: node ID of the AllocationLevelControl node +min_level: The minimum target level of the connected basin(s) +max_level: The maximum target level of the connected basin(s) +priority: If in a shortage state, the priority of the demand of the connected basin(s) +""" struct AllocationLevelControl node_id::Vector{NodeID} min_level::Vector{LinearInterpolation} From aaed23cdcb1976e5ff2ba6bcc9a26de41f7baa47 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 12 Feb 2024 15:03:04 +0100 Subject: [PATCH 23/37] Update documentation --- docs/core/allocation.qmd | 71 ++++++++++++++++++++++++++++++---------- 1 file changed, 54 insertions(+), 17 deletions(-) diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index 98c1916a8..71cb4f25a 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -56,11 +56,17 @@ $$ We consider fluxes into the basin to be positive and out of the basin to be negative. -Secondly, there is the available water in each basin above the minimum level $l_{\min,i}$ corresponding to a minimum storage $s_{\min,i}$: +Secondly, there is either a supply or demand from the storage in the basin. Given a minimum elevel $\ell_{\min, i}$ and a maximum level $\ell_{\max, i}$ which correspond to a minimum storage $s_{\min, i}$ and maximum storage $s_{\max, i}$ respectively, we get a flow supply of $$ - u_i(t)-s_{\min,i}, \quad \forall i \in B_S. + F^{\text{basin out}}_{\max, i} = \max\left(0.0, \frac{u_i(t)-s_{\max,i}}{\Delta t_{\text{alloc}}} + \phi_i(t)\right) $$ -Note that this value can be negative, which we interpret as a demand from the basin. + +and a demand of +$$ + d^p_i = \max\left(0.0, \frac{s_{\min,i} - u_i(t)}{\Delta t_{\text{alloc}}} - \phi_i(t)\right), +$$ + +for all $i \in B_S$. Here $\Delta t_{\text{alloc}}$ is the simulated time between two consecutive allocation solves. Note that the basin demand has only a single priority, so for other priorities this demand is $0$. ### Constraining factors @@ -94,47 +100,50 @@ for the set of in-neighbors and out-neighbors of a node in the allocation graph Each edge in the allocation graph has an associated capacity. These capacities are collected in the sparse capacity matrix $C_S \in \overline{\mathbb{R}}_{\ge 0}^{n'\times n'}$ where $n' = \#V'_S$ is the number of nodes in the allocation graph. The capacities can be infinite, if there is nothing in the model constraining the capacity of the edge. -The capacities are determined in 4 different ways: +The capacities are determined in different ways: - If an edge does not exist in the allocation graph, i.e. $(i,j) \notin E_S$ for certain $1 \le i,j\le n'$, then $(C_S)_{i,j} = 0$; - The capacity of the edge $e \in E_S$ is given by the smallest `max_flow_rate` of the nodes along the equivalent edges in the subnetwork. If there are no nodes with a `max_flow_rate`, the edge capacity is infinite; - If the edge is a source, the capacity of the edge is given by the flow rate of that source. +There are also capacities $C^B_S \in \mathbb{R}^b$ where $b = \# B_S$ is the number of basins, for the flow supplied by basins. + # The optimization problem The optimization problem for a subnetwork is a linear optimization problem consisting of an objective function with associated constraints on a set of variables, all of which are introduced below. ## The optimization variables -There are 2 types of variable whose value has to be determined to solve the allocation problem: +There are several types of variable whose value has to be determined to solve the allocation problem: - The flows $F \in \mathbb{R}_{\ge 0}^{n'\times n'}$ over the edges in the allocation graph; -- The flows $F^\text{basin}_{i}$ for all $i \in B_S$ supplied (positive) or consumed (negative) by the basins. +- The flows $F^\text{basin out}_{i}, F^\text{basin in}_{i} \geq 0$ for all $i \in B_S$ supplied and consumed by the basins respectively. +## The optimization objective -:::{.callout-note} -Currently the basin allocations are not taken into account in the implementation. -::: +The goal of allocation is to get the flow to nodes with demands as close as possible to these demands. To achieve this, a sum error of terms is minimized. The form of these error terms is determined by the choice of one of the supported objective function types, which are discussed further below. -## The optimization objective +$$ + \min E_{\text{user}} + E_{\text{basin}} +$$ -The goal of allocation is to get the flow to the users as close as possible to their demand. To achieve this, the following objectives are supported: +### User demands - `quadratic_absolute`: $$ - \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left( F_{ij} - d_j^p(t)\right)^2 + E_{\text{user}} = \sum_{(i,j)\in E_S\;:\; i\in U_S} \left( F_{ij} - d_j^p(t)\right)^2 $$ - `quadratic_relative`: $$ - \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left( 1 - \frac{F_{ij}}{d_j^p(t)}\right)^2 + E_{\text{user}} = \sum_{(i,j)\in E_S\;:\; i\in U_S} \left( 1 - \frac{F_{ij}}{d_j^p(t)}\right)^2 $$ - `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 + E_{\text{user}} = \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| + c \sum_{e \in E_S} F_e + E_{\text{user}} = \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 $$ :::{.callout-note} @@ -158,17 +167,45 @@ The absolute value applied here is not supported in a linear programming context In the future new optimization objectives will be introduced, for demands of basins and priorities over sources. These will be used in combination with the above, in the form of goal programming. ::: +### Basin demands + +- `quadratic_absolute`: +$$ + E_{\text{basin}} = \sum_{i \in B_S} \left( F_i^\text{basin in} - d_i^p(t)\right)^2 +$$ +- `quadratic_relative`: +$$ + E_{\text{basin}} = \sum_{i \in B_S} \left( 1 - \frac{F_i^\text{basin in}}{d_i^p(t)}\right)^2 +$$ +- `linear_absolute` (default): +$$ + E_{\text{basin}} = \sum_{i \in B_S} \left| F_i^\text{basin in} - d_i^p(t)\right| +$$ +- `linear_relative`: +$$ + E_{\text{basin}} = \sum_{i \in B_S} \left|1 - \frac{F_i^\text{basin in}}{d_i^p(t)}\right| +$$ + ## The optimization constraints - Flow conservation: For the basins in the allocation graph we have that $$ - \sum_{j=1}^{n'} F_{kj} \le F^\text{basin}_k + \sum_{i=1}^{n'} F_{ik}, \quad \forall k \in B_S . + F^\text{basin in}_k + \sum_{j=1}^{n'} F_{kj} \le F^\text{basin out}_k + \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. + +:::{.callout-note} +In @eq-flowconservationconstraint, the placement of the basin flows might seem counter-intuitive. Think of the basin storage as a separate node connected to the basin node. +::: + - 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. $$ {#eq-capacityconstraint} -By the definition of $C_S$ this also includes the source flows. +By the definition of $C_S$ this also includes the source flows. The same holds for the basin outflows: + +$$ + F^{\text{basin out}}_{i} \le F^{\text{basin out}}_{\max, i}, \quad \forall i \in B_S. +$$ :::{.callout-note} When performing subnetwork demand collection, these capacities are set to $\infty$ for edges which connect the main network to a subnetwork. From b0cfd198d0c9620a3091218ee4d40ee4aa9580c3 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 12 Feb 2024 15:09:33 +0100 Subject: [PATCH 24/37] Align implementation with documentation --- core/src/allocation_init.jl | 4 ++-- core/src/allocation_optim.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index e7c6c47df..d8a2d2b2a 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -575,11 +575,11 @@ function add_constraints_flow_conservation!( sum([ F[(node_id, outneighbor_id)] for outneighbor_id in outflow_ids_allocation(graph, node_id) - ]) + get_basin_outflow(problem, node_id) <= + ]) + get_basin_inflow(problem, node_id) <= sum([ F[(inneighbor_id, node_id)] for inneighbor_id in inflow_ids_allocation(graph, node_id) - ]) + get_basin_inflow(problem, node_id), + ]) + get_basin_outflow(problem, node_id), base_name = "flow_conservation", ) return nothing diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 37d4d4426..6b04ade03 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -448,7 +448,7 @@ function get_basin_demand( else level_min = allocation_level_control.min_level[allocation_level_control_idx](t) storage_min = get_storage_from_level(p.basin, basin_idx, level_min) - return max(0.0, (storage_basin - storage_min) / Δt_allocation - influx) + return max(0.0, (storage_min - storage_basin) / Δt_allocation - influx) end end From 02944d08edba22950dc93c66be325d9070776d42 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 13 Feb 2024 13:40:39 +0100 Subject: [PATCH 25/37] Rename AllocationLevelControl -> AllocationTarget --- core/src/allocation_optim.jl | 34 +++++++-------- core/src/parameter.jl | 6 +-- core/src/read.jl | 16 +++---- core/src/schema.jl | 8 ++-- core/src/util.jl | 13 ++---- core/src/validation.jl | 6 +-- .../AllocationLevelControlStatic.schema.json | 37 ---------------- .../AllocationLevelControlTime.schema.json | 42 ------------------- .../allocation_level_control.schema.json | 35 ---------------- python/ribasim/ribasim/__init__.py | 4 +- python/ribasim/ribasim/config.py | 14 +++---- python/ribasim/ribasim/geometry/node.py | 4 +- python/ribasim/ribasim/model.py | 6 +-- python/ribasim/ribasim/schemas.py | 4 +- .../ribasim_testmodels/__init__.py | 4 +- .../ribasim_testmodels/allocation.py | 10 ++--- .../ribasim_testmodels/basic.py | 4 +- ribasim_qgis/core/nodes.py | 12 +++--- 18 files changed, 69 insertions(+), 190 deletions(-) delete mode 100644 docs/schema/AllocationLevelControlStatic.schema.json delete mode 100644 docs/schema/AllocationLevelControlTime.schema.json delete mode 100644 docs/schema/allocation_level_control.schema.json diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 6b04ade03..a68bfb09f 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -375,7 +375,7 @@ Get several variables associated with a basin: - Its current storage - The allocation update interval - The influx (sum of instantaneous vertical fluxes of the basin) -- The index of the connected allocation_level_control node (0 if such a +- The index of the connected allocation_target node (0 if such a node does not exist) - The index of the basin """ @@ -385,27 +385,27 @@ function get_basin_data( u::ComponentVector, node_id::NodeID, ) - (; graph, basin, allocation_level_control) = p + (; graph, basin, allocation_target) = p (; Δt_allocation) = allocation_model influx = get_flow(graph, node_id, 0.0) _, basin_idx = id_index(basin.node_id, node_id) storage_basin = u.storage[basin_idx] control_inneighbors = inneighbor_labels_type(graph, node_id, EdgeType.control) if isempty(control_inneighbors) - allocation_level_control_idx = 0 + allocation_target_idx = 0 else - allocation_level_control_node_id = first(control_inneighbors) - allocation_level_control_idx = - findsorted(allocation_level_control.node_id, allocation_level_control_node_id) + allocation_target_node_id = first(control_inneighbors) + allocation_target_idx = + findsorted(allocation_target.node_id, allocation_target_node_id) end - return storage_basin, Δt_allocation, influx, allocation_level_control_idx, basin_idx + return storage_basin, Δt_allocation, influx, allocation_target_idx, basin_idx end """ Get the capacity of the basin, i.e. the maximum flow that can be abstracted from the basin if it is in a state of surplus storage (0 if no reference levels are provided by -a allocation_level_control node). +a allocation_target node). Storages are converted to flows by dividing by the allocation timestep. """ function get_basin_capacity( @@ -415,13 +415,13 @@ function get_basin_capacity( t::Float64, node_id::NodeID, )::Float64 - (; allocation_level_control) = p - storage_basin, Δt_allocation, influx, allocation_level_control_idx, basin_idx = + (; allocation_target) = p + storage_basin, Δt_allocation, influx, allocation_target_idx, basin_idx = get_basin_data(allocation_model, p, u, node_id) - if iszero(allocation_level_control_idx) + if iszero(allocation_target_idx) return 0.0 else - level_max = allocation_level_control.max_level[allocation_level_control_idx](t) + level_max = allocation_target.max_level[allocation_target_idx](t) storage_max = get_storage_from_level(p.basin, basin_idx, level_max) return max(0.0, (storage_basin - storage_max) / Δt_allocation + influx) end @@ -430,7 +430,7 @@ end """ Get the demand of the basin, i.e. how large a flow the basin needs to get to its minimum target level (0 if no -reference levels are provided by a allocation_level_control node). +reference levels are provided by a allocation_target node). Storages are converted to flows by dividing by the allocation timestep. """ function get_basin_demand( @@ -440,13 +440,13 @@ function get_basin_demand( t::Float64, node_id::NodeID, )::Float64 - (; allocation_level_control) = p - storage_basin, Δt_allocation, influx, allocation_level_control_idx, basin_idx = + (; allocation_target) = p + storage_basin, Δt_allocation, influx, allocation_target_idx, basin_idx = get_basin_data(allocation_model, p, u, node_id) - if iszero(allocation_level_control_idx) + if iszero(allocation_target_idx) return 0.0 else - level_min = allocation_level_control.min_level[allocation_level_control_idx](t) + level_min = allocation_target.min_level[allocation_target_idx](t) storage_min = get_storage_from_level(p.basin, basin_idx, level_min) return max(0.0, (storage_min - storage_basin) / Δt_allocation - influx) end diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 613daa579..b52c51b1e 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -499,12 +499,12 @@ struct User <: AbstractParameterNode end """ -node_id: node ID of the AllocationLevelControl node +node_id: node ID of the AllocationTarget node min_level: The minimum target level of the connected basin(s) max_level: The maximum target level of the connected basin(s) priority: If in a shortage state, the priority of the demand of the connected basin(s) """ -struct AllocationLevelControl +struct AllocationTarget node_id::Vector{NodeID} min_level::Vector{LinearInterpolation} max_level::Vector{LinearInterpolation} @@ -553,6 +553,6 @@ struct Parameters{T, C1, C2} discrete_control::DiscreteControl pid_control::PidControl{T} user::User - allocation_level_control::AllocationLevelControl + allocation_target::AllocationTarget subgrid::Subgrid end diff --git a/core/src/read.jl b/core/src/read.jl index f5c469d30..73cf9c344 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -764,24 +764,24 @@ function User(db::DB, config::Config)::User ) end -function AllocationLevelControl(db::DB, config::Config)::AllocationLevelControl - static = load_structvector(db, config, AllocationLevelControlStaticV1) - time = load_structvector(db, config, AllocationLevelControlTimeV1) +function AllocationTarget(db::DB, config::Config)::AllocationTarget + static = load_structvector(db, config, AllocationTargetStaticV1) + time = load_structvector(db, config, AllocationTargetTimeV1) parsed_parameters, valid = parse_static_and_time( db, config, - "AllocationLevelControl"; + "AllocationTarget"; static, time, time_interpolatables = [:min_level, :max_level], ) if !valid - error("Errors occurred when parsing AllocationLevelControl data.") + error("Errors occurred when parsing AllocationTarget data.") end - return AllocationLevelControl( + return AllocationTarget( NodeID.(parsed_parameters.node_id), parsed_parameters.min_level, parsed_parameters.max_level, @@ -874,7 +874,7 @@ function Parameters(db::DB, config::Config)::Parameters discrete_control = DiscreteControl(db, config) pid_control = PidControl(db, config, chunk_sizes) user = User(db, config) - allocation_level_control = AllocationLevelControl(db, config) + allocation_target = AllocationTarget(db, config) basin = Basin(db, config, chunk_sizes) subgrid_level = Subgrid(db, config, basin) @@ -908,7 +908,7 @@ function Parameters(db::DB, config::Config)::Parameters discrete_control, pid_control, user, - allocation_level_control, + allocation_target, subgrid_level, ) diff --git a/core/src/schema.jl b/core/src/schema.jl index c4f7fed05..0fbcf91d1 100644 --- a/core/src/schema.jl +++ b/core/src/schema.jl @@ -25,8 +25,8 @@ @schema "ribasim.outlet.static" OutletStatic @schema "ribasim.user.static" UserStatic @schema "ribasim.user.time" UserTime -@schema "ribasim.allocationlevelcontrol.static" AllocationLevelControlStatic -@schema "ribasim.allocationlevelcontrol.time" AllocationLevelControlTime +@schema "ribasim.allocationtarget.static" AllocationTargetStatic +@schema "ribasim.allocationtarget.time" AllocationTargetTime const delimiter = " / " tablename(sv::Type{SchemaVersion{T, N}}) where {T, N} = tablename(sv()) @@ -253,14 +253,14 @@ end priority::Int end -@version AllocationLevelControlStaticV1 begin +@version AllocationTargetStaticV1 begin node_id::Int min_level::Float64 max_level::Float64 priority::Int end -@version AllocationLevelControlTimeV1 begin +@version AllocationTargetTimeV1 begin node_id::Int time::DateTime min_level::Float64 diff --git a/core/src/util.jl b/core/src/util.jl index 69bc75b5d..7c0b3c5fd 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -651,24 +651,19 @@ function get_all_priorities(db::DB, config::Config)::Vector{Int} priorities = Set{Int}() # TODO: Is there a way to automatically grab all tables with a priority column? - for type in [ - UserStaticV1, - UserTimeV1, - AllocationLevelControlStaticV1, - AllocationLevelControlTimeV1, - ] + for type in [UserStaticV1, UserTimeV1, AllocationTargetStaticV1, AllocationTargetTimeV1] union!(priorities, load_structvector(db, config, type).priority) end return sort(unique(priorities)) end function get_basin_priority(p::Parameters, node_id::NodeID)::Int - (; graph, allocation_level_control) = p + (; graph, allocation_target) = p inneighbors_control = inneighbor_labels_type(graph, node_id, EdgeType.control) if isempty(inneighbors_control) return 0 else - idx = findsorted(allocation_level_control.node_id, first(inneighbors_control)) - return allocation_level_control.priority[idx] + idx = findsorted(allocation_target.node_id, first(inneighbors_control)) + return allocation_target.priority[idx] end end diff --git a/core/src/validation.jl b/core/src/validation.jl index 5e7607d52..bd7a409e1 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -3,7 +3,7 @@ neighbortypes(nodetype::Symbol) = neighbortypes(Val(nodetype)) neighbortypes(::Val{:pump}) = Set((:basin, :fractional_flow, :terminal, :level_boundary)) neighbortypes(::Val{:outlet}) = Set((:basin, :fractional_flow, :terminal, :level_boundary)) neighbortypes(::Val{:user}) = Set((:basin, :fractional_flow, :terminal, :level_boundary)) -neighbortypes(::Val{:allocation_level_control}) = Set((:basin,)) +neighbortypes(::Val{:allocation_target}) = Set((:basin,)) neighbortypes(::Val{:basin}) = Set(( :linear_resistance, :tabulated_rating_curve, @@ -58,7 +58,7 @@ n_neighbor_bounds_flow(::Val{:Terminal}) = n_neighbor_bounds(1, typemax(Int), 0, n_neighbor_bounds_flow(::Val{:PidControl}) = n_neighbor_bounds(0, 0, 0, 0) n_neighbor_bounds_flow(::Val{:DiscreteControl}) = n_neighbor_bounds(0, 0, 0, 0) n_neighbor_bounds_flow(::Val{:User}) = n_neighbor_bounds(1, 1, 1, 1) -n_neighbor_bounds_flow(::Val{:AllocationLevelControl}) = n_neighbor_bounds(0, 0, 0, 0) +n_neighbor_bounds_flow(::Val{:AllocationTarget}) = n_neighbor_bounds(0, 0, 0, 0) n_neighbor_bounds_flow(nodetype) = error("'n_neighbor_bounds_flow' not defined for $nodetype.") @@ -77,7 +77,7 @@ n_neighbor_bounds_control(::Val{:PidControl}) = n_neighbor_bounds(0, 1, 1, 1) n_neighbor_bounds_control(::Val{:DiscreteControl}) = n_neighbor_bounds(0, 0, 1, typemax(Int)) n_neighbor_bounds_control(::Val{:User}) = n_neighbor_bounds(0, 0, 0, 0) -n_neighbor_bounds_control(::Val{:AllocationLevelControl}) = +n_neighbor_bounds_control(::Val{:AllocationTarget}) = n_neighbor_bounds(0, 0, 1, typemax(Int)) n_neighbor_bounds_control(nodetype) = error("'n_neighbor_bounds_control' not defined for $nodetype.") diff --git a/docs/schema/AllocationLevelControlStatic.schema.json b/docs/schema/AllocationLevelControlStatic.schema.json deleted file mode 100644 index 18dcfa84c..000000000 --- a/docs/schema/AllocationLevelControlStatic.schema.json +++ /dev/null @@ -1,37 +0,0 @@ -{ - "$schema": "https://json-schema.org/draft/2020-12/schema", - "$id": "https://deltares.github.io/Ribasim/schema/AllocationLevelControlStatic.schema.json", - "title": "AllocationLevelControlStatic", - "description": "A AllocationLevelControlStatic object based on Ribasim.AllocationLevelControlStaticV1", - "type": "object", - "properties": { - "node_id": { - "format": "default", - "type": "integer" - }, - "min_level": { - "format": "double", - "type": "number" - }, - "max_level": { - "format": "double", - "type": "number" - }, - "priority": { - "format": "default", - "type": "integer" - }, - "remarks": { - "description": "a hack for pandera", - "type": "string", - "format": "default", - "default": "" - } - }, - "required": [ - "node_id", - "min_level", - "max_level", - "priority" - ] -} diff --git a/docs/schema/AllocationLevelControlTime.schema.json b/docs/schema/AllocationLevelControlTime.schema.json deleted file mode 100644 index 962683800..000000000 --- a/docs/schema/AllocationLevelControlTime.schema.json +++ /dev/null @@ -1,42 +0,0 @@ -{ - "$schema": "https://json-schema.org/draft/2020-12/schema", - "$id": "https://deltares.github.io/Ribasim/schema/AllocationLevelControlTime.schema.json", - "title": "AllocationLevelControlTime", - "description": "A AllocationLevelControlTime object based on Ribasim.AllocationLevelControlTimeV1", - "type": "object", - "properties": { - "node_id": { - "format": "default", - "type": "integer" - }, - "time": { - "format": "date-time", - "type": "string" - }, - "min_level": { - "format": "double", - "type": "number" - }, - "max_level": { - "format": "double", - "type": "number" - }, - "priority": { - "format": "default", - "type": "integer" - }, - "remarks": { - "description": "a hack for pandera", - "type": "string", - "format": "default", - "default": "" - } - }, - "required": [ - "node_id", - "time", - "min_level", - "max_level", - "priority" - ] -} diff --git a/docs/schema/allocation_level_control.schema.json b/docs/schema/allocation_level_control.schema.json deleted file mode 100644 index 9edae20da..000000000 --- a/docs/schema/allocation_level_control.schema.json +++ /dev/null @@ -1,35 +0,0 @@ -{ - "$schema": "https://json-schema.org/draft/2020-12/schema", - "$id": "https://deltares.github.io/Ribasim/schema/allocation_level_control.schema.json", - "title": "allocation_level_control", - "description": "A allocation_level_control object based on Ribasim.config.allocation_level_control", - "type": "object", - "properties": { - "static": { - "format": "default", - "anyOf": [ - { - "type": "null" - }, - { - "type": "string" - } - ], - "default": null - }, - "time": { - "format": "default", - "anyOf": [ - { - "type": "null" - }, - { - "type": "string" - } - ], - "default": null - } - }, - "required": [ - ] -} diff --git a/python/ribasim/ribasim/__init__.py b/python/ribasim/ribasim/__init__.py index f94ae65d3..85322bdce 100644 --- a/python/ribasim/ribasim/__init__.py +++ b/python/ribasim/ribasim/__init__.py @@ -4,7 +4,7 @@ from ribasim import utils from ribasim.config import ( Allocation, - AllocationLevelControl, + AllocationTarget, Basin, Compression, DiscreteControl, @@ -30,7 +30,7 @@ __all__ = [ "Allocation", - "AllocationLevelControl", + "AllocationTarget", "Basin", "DiscreteControl", "Compression", diff --git a/python/ribasim/ribasim/config.py b/python/ribasim/ribasim/config.py index 1456ef165..f53d21816 100644 --- a/python/ribasim/ribasim/config.py +++ b/python/ribasim/ribasim/config.py @@ -6,8 +6,8 @@ # These schemas are autogenerated from ribasim.schemas import ( - AllocationLevelControlStaticSchema, - AllocationLevelControlTimeSchema, + AllocationTargetStaticSchema, + AllocationTargetTimeSchema, BasinProfileSchema, BasinStateSchema, BasinStaticSchema, @@ -137,13 +137,13 @@ class User(NodeModel): ) -class AllocationLevelControl(NodeModel): - static: TableModel[AllocationLevelControlStaticSchema] = Field( - default_factory=TableModel[AllocationLevelControlStaticSchema], +class AllocationTarget(NodeModel): + static: TableModel[AllocationTargetStaticSchema] = Field( + default_factory=TableModel[AllocationTargetStaticSchema], json_schema_extra={"sort_keys": ["node_id", "priority"]}, ) - time: TableModel[AllocationLevelControlTimeSchema] = Field( - default_factory=TableModel[AllocationLevelControlTimeSchema], + time: TableModel[AllocationTargetTimeSchema] = Field( + default_factory=TableModel[AllocationTargetTimeSchema], json_schema_extra={"sort_keys": ["node_id", "priority", "time"]}, ) diff --git a/python/ribasim/ribasim/geometry/node.py b/python/ribasim/ribasim/geometry/node.py index 2fc31ce66..6575537e2 100644 --- a/python/ribasim/ribasim/geometry/node.py +++ b/python/ribasim/ribasim/geometry/node.py @@ -207,7 +207,7 @@ def plot(self, ax=None, zorder=None) -> Any: "DiscreteControl": "*", "PidControl": "x", "User": "s", - "AllocationLevelControl": "o", + "AllocationTarget": "o", "": "o", } @@ -225,7 +225,7 @@ def plot(self, ax=None, zorder=None) -> Any: "DiscreteControl": "k", "PidControl": "k", "User": "g", - "AllocationLevelControl": "k", + "AllocationTarget": "k", "": "k", } assert self.df is not None diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py index 21c50f951..28ff73f84 100644 --- a/python/ribasim/ribasim/model.py +++ b/python/ribasim/ribasim/model.py @@ -18,7 +18,7 @@ from ribasim.config import ( Allocation, - AllocationLevelControl, + AllocationTarget, Basin, DiscreteControl, FlowBoundary, @@ -168,9 +168,7 @@ class Model(FileModel): logging: Logging = Logging() allocation: Allocation = Field(default_factory=Allocation) - allocation_level_control: AllocationLevelControl = Field( - default_factory=AllocationLevelControl - ) + allocation_target: AllocationTarget = Field(default_factory=AllocationTarget) basin: Basin = Field(default_factory=Basin) fractional_flow: FractionalFlow = Field(default_factory=FractionalFlow) level_boundary: LevelBoundary = Field(default_factory=LevelBoundary) diff --git a/python/ribasim/ribasim/schemas.py b/python/ribasim/ribasim/schemas.py index 02565d217..1510ca320 100644 --- a/python/ribasim/ribasim/schemas.py +++ b/python/ribasim/ribasim/schemas.py @@ -11,14 +11,14 @@ class Config: coerce = True -class AllocationLevelControlStaticSchema(_BaseSchema): +class AllocationTargetStaticSchema(_BaseSchema): node_id: Series[int] = pa.Field(nullable=False) min_level: Series[float] = pa.Field(nullable=False) max_level: Series[float] = pa.Field(nullable=False) priority: Series[int] = pa.Field(nullable=False) -class AllocationLevelControlTimeSchema(_BaseSchema): +class AllocationTargetTimeSchema(_BaseSchema): node_id: Series[int] = pa.Field(nullable=False) time: Series[Timestamp] = pa.Field(nullable=False) min_level: Series[float] = pa.Field(nullable=False) diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index 42adc6c6f..107d6c3a5 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -7,7 +7,7 @@ import ribasim_testmodels from ribasim_testmodels.allocation import ( allocation_example_model, - allocation_level_control_model, + allocation_target_model, fractional_flow_subnetwork_model, looped_subnetwork_model, main_network_with_subnetworks_model, @@ -55,7 +55,7 @@ __all__ = [ "allocation_example_model", - "allocation_level_control_model", + "allocation_target_model", "backwater_model", "basic_arrow_model", "basic_model", diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 9c9ccfcc3..4d059c7a7 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -1739,20 +1739,20 @@ def main_network_with_subnetworks_model(): return model -def allocation_level_control_model(): +def allocation_target_model(): # Set up the nodes: xy = np.array( [ (0.0, 0.0), # 1: FlowBoundary (1.0, 0.0), # 2: Basin (2.0, 0.0), # 3: User - (1.0, -1.0), # 4: AllocationLevelControl + (1.0, -1.0), # 4: AllocationTarget (2.0, -1.0), # 5: Basin ] ) node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1]) - node_type = ["FlowBoundary", "Basin", "User", "AllocationLevelControl", "Basin"] + node_type = ["FlowBoundary", "Basin", "User", "AllocationTarget", "Basin"] # Make sure the feature id starts at 1: explicitly give an index. node = ribasim.Node( @@ -1810,7 +1810,7 @@ def allocation_level_control_model(): ) # Setup allocation level control - allocation_level_control = ribasim.AllocationLevelControl( + allocation_target = ribasim.AllocationTarget( static=pd.DataFrame( data={"node_id": [4], "priority": 1, "min_level": 1.0, "max_level": 1.0} ) @@ -1836,7 +1836,7 @@ def allocation_level_control_model(): network=ribasim.Network(node=node, edge=edge), basin=basin, flow_boundary=flow_boundary, - allocation_level_control=allocation_level_control, + allocation_target=allocation_target, user=user, allocation=allocation, starttime="2020-01-01 00:00:00", diff --git a/python/ribasim_testmodels/ribasim_testmodels/basic.py b/python/ribasim_testmodels/ribasim_testmodels/basic.py index 179d03a65..27011b1c8 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/basic.py +++ b/python/ribasim_testmodels/ribasim_testmodels/basic.py @@ -276,8 +276,8 @@ def basic_transient_model() -> ribasim.Model: "level": 1.4, } ) - model.basin.time = forcing # type: ignore # TODO: Fix implicit typing from pydantic. See TableModel.check_dataframe - model.basin.state = state # type: ignore # TODO: Fix implicit typing from pydantic. See TableModel.check_dataframe + model.basin.time = forcing + model.basin.state = state return model diff --git a/ribasim_qgis/core/nodes.py b/ribasim_qgis/core/nodes.py index 1c034e857..62b28b37e 100644 --- a/ribasim_qgis/core/nodes.py +++ b/ribasim_qgis/core/nodes.py @@ -215,9 +215,9 @@ def renderer(self) -> QgsCategorizedSymbolRenderer: "DiscreteControl": (QColor("black"), "DiscreteControl", shape.Star), "PidControl": (QColor("black"), "PidControl", shape.Cross2), "User": (QColor("green"), "User", shape.Square), - "AllocationLevelControl": ( + "AllocationTarget": ( QColor("black"), - "AllocationLevelControl", + "AllocationTarget", shape.Circle, ), # All other nodes, or incomplete input @@ -785,10 +785,10 @@ def attributes(cls) -> list[QgsField]: ] -class AllocationLevelControl(Input): +class AllocationTarget(Input): @classmethod def input_type(cls) -> str: - return "AllocationLevelControl / static" + return "AllocationTarget / static" @classmethod def geometry_type(cls) -> str: @@ -803,10 +803,10 @@ def attributes(cls) -> list[QgsField]: ] -class AllocationLevelControlTime(Input): +class AllocationTargetTime(Input): @classmethod def input_type(cls) -> str: - return "AllocationLevelControl / time" + return "AllocationTarget / time" @classmethod def geometry_type(cls) -> str: From 724caa0d8400331a0d69f3184e65eb45efc387ef Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 13 Feb 2024 13:58:29 +0100 Subject: [PATCH 26/37] Merge fix --- core/test/allocation_test.jl | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index c4494ddb4..52b7cedce 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -102,8 +102,8 @@ end F = problem[:F] F_abs_user = problem[:F_abs_user] - @test objective.terms[F_abs[NodeID(:User, 5)]] == 1.0 - @test objective.terms[F_abs[NodeID(:User, 6)]] == 1.0 + @test objective.terms[F_abs_user[NodeID(:User, 5)]] == 1.0 + @test objective.terms[F_abs_user[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 @@ -119,8 +119,8 @@ end F = problem[:F] F_abs_user = problem[:F_abs_user] - @test objective.terms[F_abs[NodeID(:User, 5)]] == 1.0 - @test objective.terms[F_abs[NodeID(:User, 6)]] == 1.0 + @test objective.terms[F_abs_user[NodeID(:User, 5)]] == 1.0 + @test objective.terms[F_abs_user[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 @@ -221,9 +221,9 @@ end # support absolute value expressions in the objective function allocation_model_main_network = Ribasim.get_allocation_model(p, 1) problem = allocation_model_main_network.problem - @test problem[:F_abs].axes[1] == NodeID.(:Pump, [11, 24, 38]) - @test problem[:abs_positive].axes[1] == NodeID.(:Pump, [11, 24, 38]) - @test problem[:abs_negative].axes[1] == NodeID.(:Pump, [11, 24, 38]) + @test problem[:F_abs_user].axes[1] == NodeID.(:Pump, [11, 24, 38]) + @test problem[:abs_positive_user].axes[1] == NodeID.(:Pump, [11, 24, 38]) + @test problem[:abs_negative_user].axes[1] == NodeID.(:Pump, [11, 24, 38]) # In each subnetwork, the connection from the main network to the subnetwork is # interpreted as a source @@ -277,14 +277,15 @@ end # Main network objective function objective = JuMP.objective_function(problem) objective_variables = keys(objective.terms) - F_abs = problem[:F_abs] - @test F_abs[NodeID(:Pump, 11)] ∈ objective_variables - @test F_abs[NodeID(:Pump, 24)] ∈ objective_variables - @test F_abs[NodeID(:Pump, 38)] ∈ objective_variables + F_abs_user = problem[:F_abs_user] + @test F_abs_user[NodeID(:Pump, 11)] ∈ objective_variables + @test F_abs_user[NodeID(:Pump, 24)] ∈ objective_variables + @test F_abs_user[NodeID(:Pump, 38)] ∈ objective_variables # Running full allocation algorithm Ribasim.set_flow!(graph, NodeID(:FlowBoundary, 1), NodeID(:Basin, 2), 4.5) - Ribasim.update_allocation!((; p, t)) + u = ComponentVector(; storage = zeros(length(p.basin.node_id))) + Ribasim.update_allocation!((; p, t, u)) @test subnetwork_allocateds[NodeID(:Basin, 2), NodeID(:Pump, 11)] ≈ [4.0, 0.49766666, 0.0] From ae49fb2c00f2c97968d777d451720ae865f36007 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 13 Feb 2024 16:00:42 +0100 Subject: [PATCH 27/37] happify mypy --- python/ribasim_testmodels/ribasim_testmodels/basic.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/ribasim_testmodels/ribasim_testmodels/basic.py b/python/ribasim_testmodels/ribasim_testmodels/basic.py index 27011b1c8..179d03a65 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/basic.py +++ b/python/ribasim_testmodels/ribasim_testmodels/basic.py @@ -276,8 +276,8 @@ def basic_transient_model() -> ribasim.Model: "level": 1.4, } ) - model.basin.time = forcing - model.basin.state = state + model.basin.time = forcing # type: ignore # TODO: Fix implicit typing from pydantic. See TableModel.check_dataframe + model.basin.state = state # type: ignore # TODO: Fix implicit typing from pydantic. See TableModel.check_dataframe return model From d1ae0f8289b02de6c0301bb56f48ae22b0ec99a9 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 13 Feb 2024 16:20:15 +0100 Subject: [PATCH 28/37] Fix initialization of AllocationTarget --- core/src/read.jl | 2 +- core/test/allocation_test.jl | 8 ++++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/core/src/read.jl b/core/src/read.jl index 56d87aebc..a6ed7e595 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -779,7 +779,7 @@ function AllocationTarget(db::DB, config::Config)::AllocationTarget end return AllocationTarget( - NodeID.(parsed_parameters.node_id), + NodeID.(NodeType.AllocationTarget, parsed_parameters.node_id), parsed_parameters.min_level, parsed_parameters.max_level, parsed_parameters.priority, diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 9d9658a88..eaf54c3e8 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -319,3 +319,11 @@ end @test user.allocated[2] ≈ [4.0, 0.0, 0.0] @test user.allocated[7] ≈ [0.001, 0.0, 0.0] end + +@testitem "Allocation level control" begin + toml_path = + normpath(@__DIR__, "../../generated_testmodels/allocation_target/ribasim.toml") + @test ispath(toml_path) + + Ribasim.run(toml_path) +end From fcee66a3cb4e3c7cdd714b27e8046c6e3ae463d8 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 14 Feb 2024 09:20:49 +0100 Subject: [PATCH 29/37] test model adjustments --- .../ribasim_testmodels/allocation.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 4d059c7a7..982ac79bb 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -1794,10 +1794,10 @@ def allocation_target_model(): static = pd.DataFrame( data={ "node_id": [2, 5], - "drainage": 0.1e-3, - "potential_evaporation": 0.2e-3, - "infiltration": 0.3e-3, - "precipitation": 0.4e-3, + "drainage": 0.0, + "potential_evaporation": 0.0, + "infiltration": 0.0, + "precipitation": 1e-6, "urban_runoff": 0.0, } ) @@ -1812,7 +1812,7 @@ def allocation_target_model(): # Setup allocation level control allocation_target = ribasim.AllocationTarget( static=pd.DataFrame( - data={"node_id": [4], "priority": 1, "min_level": 1.0, "max_level": 1.0} + data={"node_id": [4], "priority": 1, "min_level": 1.0, "max_level": 1.5} ) ) @@ -1822,7 +1822,7 @@ def allocation_target_model(): data={ "node_id": [3], "priority": [2], - "demand": [1e-3], + "demand": [1.5e-3], "return_factor": [0.2], "min_level": [0.2], } From 9a91ce1d731b96b45e09d2cc748f470ac1014fa7 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 14 Feb 2024 11:02:52 +0100 Subject: [PATCH 30/37] Add tests --- core/test/allocation_test.jl | 36 ++++++++++++++++++- .../ribasim_testmodels/allocation.py | 4 +-- 2 files changed, 37 insertions(+), 3 deletions(-) diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index eaf54c3e8..47fec1842 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -324,6 +324,40 @@ end toml_path = normpath(@__DIR__, "../../generated_testmodels/allocation_target/ribasim.toml") @test ispath(toml_path) + model = Ribasim.run(toml_path) - Ribasim.run(toml_path) + storage = Ribasim.get_storages_and_levels(model).storage[1, :] + t = Ribasim.timesteps(model) + + p = model.integrator.p + (; user, graph, allocation, basin, allocation_target) = p + + d = user.demand_itp[1][2](0) + ϕ = Ribasim.get_flow(graph, Ribasim.NodeID(Ribasim.NodeType.Basin, 2), 0) + q = Ribasim.get_flow( + graph, + Ribasim.NodeID(Ribasim.NodeType.FlowBoundary, 1), + Ribasim.NodeID(Ribasim.NodeType.Basin, 2), + 0, + ) + A = basin.area[1][1] + l_max = allocation_target.max_level[1](0) + Δt_allocation = allocation.allocation_models[1].Δt_allocation + + # Until the first allocation solve, the user abstracts fully + pre_allocation = t .<= Δt_allocation + u_pre_allocation(τ) = storage[1] + (q + ϕ - d) * τ + @test storage[pre_allocation] ≈ u_pre_allocation.(t[pre_allocation]) atol = 4e-3 + + # Until the basin is at its maximum level, the user does not abstract + basin_filling = @. ~pre_allocation && (storage <= A * l_max) + fill_start_idx = findlast(pre_allocation) + u_filling(τ) = storage[fill_start_idx] + (q + ϕ) * (τ - t[fill_start_idx]) + @test storage[basin_filling] ≈ u_filling.(t[basin_filling]) atol = 1 + + # After the basin has reached its maximum level, the user abstracts fully again + after_filling = @. ~pre_allocation && ~basin_filling + fill_stop_idx = findfirst(after_filling) + u_after_filling(τ) = storage[fill_stop_idx] + (q + ϕ - d) * (τ - t[fill_stop_idx]) + @test storage[after_filling] ≈ u_after_filling.(t[after_filling]) atol = 4e-1 end diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 982ac79bb..a7fc3496c 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -1830,7 +1830,7 @@ def allocation_target_model(): ) # Setup allocation - allocation = ribasim.Allocation(use_allocation=True, timestep=86400) + allocation = ribasim.Allocation(use_allocation=True, timestep=1e5) model = ribasim.Model( network=ribasim.Network(node=node, edge=edge), @@ -1840,7 +1840,7 @@ def allocation_target_model(): user=user, allocation=allocation, starttime="2020-01-01 00:00:00", - endtime="2020-03-01 00:00:00", + endtime="2020-01-21 00:00:00", ) return model From 1f0876b5690ed3597a111f6880e79b722924ff5e Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 14 Feb 2024 13:00:08 +0100 Subject: [PATCH 31/37] Small fixes --- core/src/allocation_init.jl | 8 +++++--- core/test/allocation_test.jl | 6 +++--- .../ribasim_testmodels/allocation.py | 20 +++++++++++++++---- 3 files changed, 24 insertions(+), 10 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 11d055cdb..3d8148ced 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -560,6 +560,8 @@ function add_constraints_flow_conservation!( )::Nothing (; graph) = p F = problem[:F] + F_basin_in = problem[:F_basin_in] + F_basin_out = problem[:F_basin_out] node_ids = graph[].node_ids[allocation_network_id] node_ids_conservation = [node_id for node_id in node_ids if node_id.type == NodeType.Basin] @@ -571,14 +573,14 @@ function add_constraints_flow_conservation!( problem[:flow_conservation] = JuMP.@constraint( problem, [node_id = node_ids_conservation], - sum([ + get_basin_inflow(problem, node_id) + sum([ F[(node_id, outneighbor_id)] for outneighbor_id in outflow_ids_allocation(graph, node_id) ]) == - sum([ + get_basin_outflow(problem, node_id) + sum([ F[(inneighbor_id, node_id)] for inneighbor_id in inflow_ids_allocation(graph, node_id) - ]) + get_basin_outflow(problem, node_id), + ]), base_name = "flow_conservation", ) return nothing diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 47fec1842..86893ef0c 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -347,17 +347,17 @@ end # Until the first allocation solve, the user abstracts fully pre_allocation = t .<= Δt_allocation u_pre_allocation(τ) = storage[1] + (q + ϕ - d) * τ - @test storage[pre_allocation] ≈ u_pre_allocation.(t[pre_allocation]) atol = 4e-3 + @test storage[pre_allocation] ≈ u_pre_allocation.(t[pre_allocation]) rtol = 1e-4 # Until the basin is at its maximum level, the user does not abstract basin_filling = @. ~pre_allocation && (storage <= A * l_max) fill_start_idx = findlast(pre_allocation) u_filling(τ) = storage[fill_start_idx] + (q + ϕ) * (τ - t[fill_start_idx]) - @test storage[basin_filling] ≈ u_filling.(t[basin_filling]) atol = 1 + @test storage[basin_filling] ≈ u_filling.(t[basin_filling]) rtol = 1e-4 # After the basin has reached its maximum level, the user abstracts fully again after_filling = @. ~pre_allocation && ~basin_filling fill_stop_idx = findfirst(after_filling) u_after_filling(τ) = storage[fill_stop_idx] + (q + ϕ - d) * (τ - t[fill_stop_idx]) - @test storage[after_filling] ≈ u_after_filling.(t[after_filling]) atol = 4e-1 + @test storage[after_filling] ≈ u_after_filling.(t[after_filling]) rtol = 1e-4 end diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index a7fc3496c..b8a2fb18b 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -1793,16 +1793,28 @@ def allocation_target_model(): ) static = pd.DataFrame( data={ - "node_id": [2, 5], + "node_id": [5], "drainage": 0.0, "potential_evaporation": 0.0, "infiltration": 0.0, - "precipitation": 1e-6, + "precipitation": 0.0, "urban_runoff": 0.0, } ) + time = pd.DataFrame( + data={ + "node_id": 2, + "time": ["2020-01-01 00:00:00", "2020-01-16 00:00:00"], + "drainage": 0.0, + "potential_evaporation": 0.0, + "infiltration": 0.0, + "precipitation": [1e-6, 0.0], + "urban_runoff": 0.0, + }, + ) + state = pd.DataFrame(data={"node_id": [2, 5], "level": 0.5}) - basin = ribasim.Basin(profile=profile, static=static, state=state) + basin = ribasim.Basin(profile=profile, static=static, time=time, state=state) # Setup flow boundary flow_boundary = ribasim.FlowBoundary( @@ -1840,7 +1852,7 @@ def allocation_target_model(): user=user, allocation=allocation, starttime="2020-01-01 00:00:00", - endtime="2020-01-21 00:00:00", + endtime="2020-02-01 00:00:00", ) return model From 82b566049da3db8928f40d8c3e9c5904c494d871 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 14 Feb 2024 14:55:36 +0100 Subject: [PATCH 32/37] Finalize tests --- core/test/allocation_test.jl | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 86893ef0c..e9ebac9b7 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -333,7 +333,7 @@ end (; user, graph, allocation, basin, allocation_target) = p d = user.demand_itp[1][2](0) - ϕ = Ribasim.get_flow(graph, Ribasim.NodeID(Ribasim.NodeType.Basin, 2), 0) + ϕ = 1e-3 # precipitation q = Ribasim.get_flow( graph, Ribasim.NodeID(Ribasim.NodeType.FlowBoundary, 1), @@ -356,8 +356,17 @@ end @test storage[basin_filling] ≈ u_filling.(t[basin_filling]) rtol = 1e-4 # After the basin has reached its maximum level, the user abstracts fully again - after_filling = @. ~pre_allocation && ~basin_filling + precipitation = eachindex(storage) .<= argmax(storage) + after_filling = @. ~pre_allocation && ~basin_filling && precipitation fill_stop_idx = findfirst(after_filling) u_after_filling(τ) = storage[fill_stop_idx] + (q + ϕ - d) * (τ - t[fill_stop_idx]) @test storage[after_filling] ≈ u_after_filling.(t[after_filling]) rtol = 1e-4 + + # After precipitation stops, the user still abstracts from the basin so the storage decreases + storage_reduction = @. ~precipitation && t <= 1.8e6 + storage_reduction_start = findfirst(storage_reduction) + u_storage_reduction(τ) = + storage[storage_reduction_start] + (q - d) * (τ - t[storage_reduction_start]) + @test storage[storage_reduction] ≈ u_storage_reduction.(t[storage_reduction]) rtol = + 1e-4 end From 8fd9d5c47f11d7dfdfa718ba772df8e8a0463558 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 14 Feb 2024 17:11:27 +0100 Subject: [PATCH 33/37] Add AllocationTarget node to usage.qmd --- docs/core/usage.qmd | 33 +++++++++++++++++++++++++++++++-- 1 file changed, 31 insertions(+), 2 deletions(-) diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index 438115565..e06645c9d 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -198,6 +198,9 @@ name it must have in the database if it is stored there. - User: sets water usage demands at certain priorities - `User / static`: demands - `User / time`: dynamic demands +- AllocationTarget: Indicates minimum and maximum and maximum target level of connected basins for allocation + - `AllocationTarget / static`: static target levels + - `AllocationTarget / time`: dynamic target levels - Terminal: Water sink without state or properties - `Terminal / static`: - (only node IDs) - DiscreteControl: Set parameters of other nodes based on model state conditions (e.g. basin level) @@ -471,7 +474,7 @@ active | Bool | - | (optional, default true) demand | Float64 | $m^3 s^{-1}$ | - return_factor | Float64 | - | between [0 - 1] min_level | Float64 | $m$ | - -priority | Int | - | sorted per node id +priority | Int | - | positive, sorted per node id ## User / time @@ -486,12 +489,38 @@ Note that a `node_id` can be either in this table or in the static one, but not column | type | unit | restriction ------------- | -------- | ------------ | ----------- node_id | Int | - | sorted -priority | Int | - | sorted per node id +priority | Int | - | positive, sorted per node id time | DateTime | - | sorted per priority per node id demand | Float64 | $m^3 s^{-1}$ | - return_factor | Float64 | - | between [0 - 1] min_level | Float64 | $m$ | (optional) +# AllocationTarget {#sec-allocation_target} + +An `AllocationTarget` 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. + +column | type | unit | restriction +------------- | ------- | ------------ | ----------- +node_id | Int | - | sorted +min_level | Float64 | $m$ | - +max_level | Float64 | $m$ | - +priority | Int | - | positive + +## AllocationTarget / time + +This table is the transient form of the `AllocationTarget` table, in which time-dependent minimum and maximum levels can be supplied. +Similar to the static version, only a single priority per `AllocationTarget` node can be provided. + +column | type | unit | restriction +------------- | ------- | ------------ | ----------- +node_id | Int | - | sorted +time | DateTime | - | sorted per node id +min_level | Float64 | $m$ | - +max_level | Float64 | $m$ | - +priority | Int | - | positive + ## Allocation results The allocation table contains a record of allocation results: when it happened, for which user node, in which allocation network, and what the demand, allocated flow and abstracted flow were. From 4e7d8e50659bc1f10295c694388795015295f35a Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 14 Feb 2024 18:24:52 +0100 Subject: [PATCH 34/37] Introduce results section in usage.qmd --- docs/core/usage.qmd | 139 +++++++++++++++++++++++++------------------- 1 file changed, 79 insertions(+), 60 deletions(-) diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index e06645c9d..bfd1c8ad8 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -345,37 +345,6 @@ Water levels beyond the last `basin_level` are linearly extrapolated. Note that the interpolation to subgrid water level is not constrained by any water balance within Ribasim. Generally, to create physically meaningful subgrid water levels, the subgrid table must be parametrized properly such that the spatially integrated water volume of the subgrid elements agrees with the total storage volume of the basin. -## Basin results - -The basin table contains results of the storage and level of each basin at every solver -timestep. The initial condition is also written to the file. - -column | type | unit --------- | -------- | ---- -time | DateTime | - -node_id | Int | - -storage | Float64 | $m^3$ -level | Float64 | $m$ - -The table is sorted by time, and per time it is sorted by `node_id`. - -## Flow results - -The flow table contains results of the flow on every edge in the model, for each solver -timestep. - -column | type | unit -------------- | ------------------- | ---- -time | DateTime | - -edge_id | Union{Int, Missing} | - -from_node_id | Int | - -to_node_id | Int | - -flow | Float64 | $m^3 s^{-1}$ - -The table is sorted by time, and per time the same `edge_id` order is used, though not sorted. -Flows that are added to the model at a node, have a missing `edge_id`, and identical `from_node_id` and `to_node_id`. -Flows out of the model always have a negative sign, and additions a positive sign. - # FractionalFlow {#sec-fractional-flow} Lets a fraction (in [0,1]) of the incoming flow trough. @@ -521,24 +490,6 @@ min_level | Float64 | $m$ | - max_level | Float64 | $m$ | - priority | Int | - | positive -## Allocation results - -The allocation table contains a record of allocation results: when it happened, for which user node, in which allocation network, and what the demand, allocated flow and abstracted flow were. - -column | type ---------------------- | ----- -time | DateTime -allocation_network_id | Int -user_node_id | Int -priority | Int -demand | Float64 -allocated | Float64 -abstracted | Float64 - -::: {.callout-note} -Currently the stored demand and abstraction rate are those at the allocation timepoint (and the abstraction rate is based on the previous allocation optimization). In the future these will be an average over the previous allocation timestep. -::: - # LevelBoundary {#sec-level-boundary} Acts like an infinitely large basin where the level does not change by flow. @@ -671,17 +622,6 @@ node_id | Int | - | sorted control_state | String | - | - truth_state | String | - | Consists of the characters "T" (true), "F" (false), "U" (upcrossing), "D" (downcrossing) and "*" (any), sorted per node_id -## DiscreteControl results - -The control table contains a record of each change of control state: when it happened, which control node was involved, to which control state it changed and based on which truth state. - -column | type ---------------- | ---------- -time | DateTime -control_node_id | Int -truth_state | String -control_state | String - # PidControl {#sec-pid-control} The PidControl node controls the level in a basin by continuously controlling the flow rate of a connected pump or outlet. See also [PID controller](https://en.wikipedia.org/wiki/PID_controller). When A PidControl node is made inactive, the node under its control retains the last flow rate value, and the error integral is reset to 0. @@ -718,3 +658,82 @@ target | Float64 | $m$ | - proportional | Float64 | $s^{-1}$ | - integral | Float64 | $s^{-2}$ | - derivative | Float64 | - | - + +# Results + +## Basin (storage) + +The basin table contains results of the storage and level of each basin at every solver +timestep. The initial condition is also written to the file. + +column | type | unit +-------- | -------- | ---- +time | DateTime | - +node_id | Int | - +storage | Float64 | $m^3$ +level | Float64 | $m$ + +The table is sorted by time, and per time it is sorted by `node_id`. + +## Flow + +The flow table contains results of the flow on every edge in the model, for each solver +timestep. + +column | type | unit +------------- | ------------------- | ---- +time | DateTime | - +edge_id | Union{Int, Missing} | - +from_node_id | Int | - +to_node_id | Int | - +flow | Float64 | $m^3 s^{-1}$ + +The table is sorted by time, and per time the same `edge_id` order is used, though not sorted. +Flows that are added to the model at a node, have a missing `edge_id`, and identical `from_node_id` and `to_node_id`. +Flows out of the model always have a negative sign, and additions a positive sign. + +## Control (DiscreteControl) + +The control table contains a record of each change of control state: when it happened, which control node was involved, to which control state it changed and based on which truth state. + +column | type +--------------- | ---------- +time | DateTime +control_node_id | Int +truth_state | String +control_state | String + +## Allocation + +The allocation table contains a record of allocation results: when it happened, for which user node, in which allocation network, and what the demand, allocated flow and abstracted flow were. + +column | type +--------------------- | ----- +time | DateTime +allocation_network_id | Int +user_node_id | Int +priority | Int +demand | Float64 +allocated | Float64 +abstracted | Float64 + +::: {.callout-note} +Currently the stored demand and abstraction rate are those at the allocation timepoint (and the abstraction rate is based on the previous allocation optimization). In the future these will be an average over the previous allocation timestep. +::: + +## Allocation flow + +The allocation flow table contains results of the optimized allocation flow on every edge in the model that is part of a subnetwork, for each time an optimization problem is solved (see also [here](allocation.qmd#the-high-level-algorithm)). +If in the model a main network and subnetwork(s) are specified, there are 2 different +types of optimization for the subnetwork: collecting its total demand per priority (for allocating flow from the main network to the subnetwork), and allocating flow within the subnetwork. The column `collect_demands` provides the distinction between these two optimization types. + +column | type +--------------------- | ----- +time | DateTime +edge_id | Int +from_node_id | Int +to_node_id | Int +allocation_network_id | Int +priority | Int +flow | Float64 +collect_demands | Bool From 46550470c4e410f64c91d636d9e5a7a6854dc63e Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 16 Feb 2024 10:55:45 +0100 Subject: [PATCH 35/37] Comments adressed --- core/src/allocation_init.jl | 6 +----- core/src/allocation_optim.jl | 7 +++++-- core/src/parameter.jl | 9 ++++----- core/src/util.jl | 15 +++++++++------ core/test/allocation_test.jl | 4 ++-- docs/core/allocation.qmd | 2 +- docs/core/usage.qmd | 14 +++++++------- ribasim_qgis/core/nodes.py | 2 +- 8 files changed, 30 insertions(+), 29 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 3d8148ced..21ed1eb5b 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -560,8 +560,6 @@ function add_constraints_flow_conservation!( )::Nothing (; graph) = p F = problem[:F] - F_basin_in = problem[:F_basin_in] - F_basin_out = problem[:F_basin_out] node_ids = graph[].node_ids[allocation_network_id] node_ids_conservation = [node_id for node_id in node_ids if node_id.type == NodeType.Basin] @@ -782,13 +780,11 @@ end """ Add the basin flow constraints to the allocation problem. The constraint indices are the basin node IDs. -NOTE: negative basin flow means flow out of the basin. Constraint: -flow basin >= - basin capacity +flow out of basin <= basin capacity """ function add_constraints_basin_flow!(problem::JuMP.Model)::Nothing - F_basin_in = problem[:F_basin_in] F_basin_out = problem[:F_basin_out] problem[:basin_outflow] = JuMP.@constraint( problem, diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 184431332..76c0eff5b 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -163,9 +163,9 @@ function set_objective_priority!( user_idx = findsorted(node_id, node_id_user) if demand_from_timeseries[user_idx] d = demand_itp[user_idx][priority_idx](t) - set_user_demand!(user, node_id_user, priority_idx, d) + set_user_demand!(p, node_id_user, priority_idx, d) else - d = get_user_demand(user, node_id_user, priority_idx) + d = get_user_demand(p, node_id_user, priority_idx) end demand_max = max(demand_max, d) add_user_term!(ex, edge_id, objective_type, d, problem) @@ -372,6 +372,7 @@ function get_basin_data( ) (; graph, basin, allocation_target) = p (; Δt_allocation) = allocation_model + @assert node_id.type == NodeType.Basin influx = get_flow(graph, node_id, 0.0) _, basin_idx = id_index(basin.node_id, node_id) storage_basin = u.storage[basin_idx] @@ -401,6 +402,7 @@ function get_basin_capacity( node_id::NodeID, )::Float64 (; allocation_target) = p + @assert node_id.type == NodeType.Basin storage_basin, Δt_allocation, influx, allocation_target_idx, basin_idx = get_basin_data(allocation_model, p, u, node_id) if iszero(allocation_target_idx) @@ -426,6 +428,7 @@ function get_basin_demand( node_id::NodeID, )::Float64 (; allocation_target) = p + @assert node_id.type == NodeType.Basin storage_basin, Δt_allocation, influx, allocation_target_idx, basin_idx = get_basin_data(allocation_model, p, u, node_id) if iszero(allocation_target_idx) diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 2f6d87bb9..c7e8ae6b4 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -63,6 +63,7 @@ allocation_network_ids: The unique sorted allocation network IDs allocation models: The allocation models for the main network and subnetworks corresponding to allocation_network_ids main_network_connections: (from_id, to_id) from the main network to the subnetwork per subnetwork +priorities: All used priority values. subnetwork_demands: The demand of an edge from the main network to a subnetwork record: A record of all flows computed by allocation optimization, eventually saved to output file @@ -457,13 +458,13 @@ struct PidControl{T} <: AbstractParameterNode end """ -demand: water flux demand of user per priority over time +demand: water flux demand of user per priority over time. + Each user has a demand for all priorities, + which is 0.0 if it is not provided explicitly. active: whether this node is active and thus demands water allocated: water flux currently allocated to user per priority return_factor: the factor in [0,1] of how much of the abstracted water is given back to the system min_level: The level of the source basin below which the user does not abstract -priorities: All used priority values. Each user has a demand for all these priorities, - which is 0.0 if it is not provided explicitly. record: Collected data of allocation optimizations for output file. """ struct User <: AbstractParameterNode @@ -475,7 +476,6 @@ struct User <: AbstractParameterNode allocated::Vector{Vector{Float64}} return_factor::Vector{Float64} min_level::Vector{Float64} - priorities::Vector{Int} record::@NamedTuple{ time::Vector{Float64}, allocation_network_id::Vector{Int}, @@ -508,7 +508,6 @@ struct User <: AbstractParameterNode allocated, return_factor, min_level, - priorities, record, ) else diff --git a/core/src/util.jl b/core/src/util.jl index be276bdb5..8b1d4ac8e 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -378,7 +378,7 @@ function get_level( storage::Union{AbstractArray, Number} = 0, )::Union{Real, Nothing} (; basin, level_boundary) = p - if node_id.type == NodeType.Basin + if node_id.type == NodeType.Basin _, i = id_index(basin.node_id, node_id) current_level = get_tmp(basin.current_level, storage) current_level[i] @@ -620,22 +620,24 @@ function is_main_network(allocation_network_id::Int)::Bool return allocation_network_id == 1 end -function get_user_demand(user::User, node_id::NodeID, priority_idx::Int)::Float64 +function get_user_demand(p::Parameters, node_id::NodeID, priority_idx::Int)::Float64 + (; user, allocation) = p (; demand) = user user_idx = findsorted(user.node_id, node_id) - n_priorities = length(user.priorities) + n_priorities = length(allocation.priorities) return demand[(user_idx - 1) * n_priorities + priority_idx] end function set_user_demand!( - user::User, + p::Parameters, node_id::NodeID, priority_idx::Int, value::Float64, )::Nothing + (; user, allocation) = p (; demand) = user user_idx = findsorted(user.node_id, node_id) - n_priorities = length(user.priorities) + n_priorities = length(allocation.priorities) demand[(user_idx - 1) * n_priorities + priority_idx] = value return nothing end @@ -652,11 +654,12 @@ end function get_basin_priority(p::Parameters, node_id::NodeID)::Int (; graph, allocation_target) = p + @assert node_id.type == NodeType.Basin inneighbors_control = inneighbor_labels_type(graph, node_id, EdgeType.control) if isempty(inneighbors_control) return 0 else - idx = findsorted(allocation_target.node_id, first(inneighbors_control)) + idx = findsorted(allocation_target.node_id, only(inneighbors_control)) return allocation_target.priority[idx] end end diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index e9ebac9b7..e99909267 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -44,9 +44,9 @@ # Test getting and setting user demands (; user) = p - Ribasim.set_user_demand!(user, NodeID(:User, 11), 2, Float64(π)) + Ribasim.set_user_demand!(p, NodeID(:User, 11), 2, Float64(π)) @test user.demand[4] ≈ π - @test Ribasim.get_user_demand(user, NodeID(:User, 11), 2) ≈ π + @test Ribasim.get_user_demand(p, NodeID(:User, 11), 2) ≈ π end @testitem "Allocation objective: quadratic absolute" skip = true begin diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index 8ce6bc580..d3303264b 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -56,7 +56,7 @@ $$ We consider fluxes into the basin to be positive and out of the basin to be negative. -Secondly, there is either a supply or demand from the storage in the basin. Given a minimum elevel $\ell_{\min, i}$ and a maximum level $\ell_{\max, i}$ which correspond to a minimum storage $s_{\min, i}$ and maximum storage $s_{\max, i}$ respectively, we get a flow supply of +Secondly, there is either a supply or demand from the storage in the basin. Given a minimum level $\ell_{\min, i}$ and a maximum level $\ell_{\max, i}$ which correspond to a minimum storage $s_{\min, i}$ and maximum storage $s_{\max, i}$ respectively, we get a flow supply of $$ F^{\text{basin out}}_{\max, i} = \max\left(0.0, \frac{u_i(t)-s_{\max,i}}{\Delta t_{\text{alloc}}} + \phi_i(t)\right) $$ diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index bfd1c8ad8..9d699d822 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -198,7 +198,7 @@ name it must have in the database if it is stored there. - User: sets water usage demands at certain priorities - `User / static`: demands - `User / time`: dynamic demands -- AllocationTarget: Indicates minimum and maximum and maximum target level of connected basins for allocation +- AllocationTarget: Indicates minimum and maximum target level of connected basins for allocation - `AllocationTarget / static`: static target levels - `AllocationTarget / time`: dynamic target levels - Terminal: Water sink without state or properties @@ -462,7 +462,7 @@ priority | Int | - | positive, sorted per node id time | DateTime | - | sorted per priority per node id demand | Float64 | $m^3 s^{-1}$ | - return_factor | Float64 | - | between [0 - 1] -min_level | Float64 | $m$ | (optional) +min_level | Float64 | $m$ | - # AllocationTarget {#sec-allocation_target} @@ -661,7 +661,7 @@ derivative | Float64 | - | - # Results -## Basin (storage) +## Basin - `basin.arrow` The basin table contains results of the storage and level of each basin at every solver timestep. The initial condition is also written to the file. @@ -675,7 +675,7 @@ level | Float64 | $m$ The table is sorted by time, and per time it is sorted by `node_id`. -## Flow +## Flow - `flow.arrow` The flow table contains results of the flow on every edge in the model, for each solver timestep. @@ -692,7 +692,7 @@ The table is sorted by time, and per time the same `edge_id` order is used, thou Flows that are added to the model at a node, have a missing `edge_id`, and identical `from_node_id` and `to_node_id`. Flows out of the model always have a negative sign, and additions a positive sign. -## Control (DiscreteControl) +## DiscreteControl - `control.arrow` The control table contains a record of each change of control state: when it happened, which control node was involved, to which control state it changed and based on which truth state. @@ -703,7 +703,7 @@ control_node_id | Int truth_state | String control_state | String -## Allocation +## Allocation - `allocation.arrow` The allocation table contains a record of allocation results: when it happened, for which user node, in which allocation network, and what the demand, allocated flow and abstracted flow were. @@ -721,7 +721,7 @@ abstracted | Float64 Currently the stored demand and abstraction rate are those at the allocation timepoint (and the abstraction rate is based on the previous allocation optimization). In the future these will be an average over the previous allocation timestep. ::: -## Allocation flow +## Allocation flow - `allocation_flow.arrow` The allocation flow table contains results of the optimized allocation flow on every edge in the model that is part of a subnetwork, for each time an optimization problem is solved (see also [here](allocation.qmd#the-high-level-algorithm)). If in the model a main network and subnetwork(s) are specified, there are 2 different diff --git a/ribasim_qgis/core/nodes.py b/ribasim_qgis/core/nodes.py index 62b28b37e..9e1ddd79f 100644 --- a/ribasim_qgis/core/nodes.py +++ b/ribasim_qgis/core/nodes.py @@ -785,7 +785,7 @@ def attributes(cls) -> list[QgsField]: ] -class AllocationTarget(Input): +class AllocationTargetStatic(Input): @classmethod def input_type(cls) -> str: return "AllocationTarget / static" From f0e70c898470658300676fb6322361248382f5a6 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 16 Feb 2024 11:20:42 +0100 Subject: [PATCH 36/37] Fix saving results --- core/src/allocation_optim.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 76c0eff5b..812875639 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -228,7 +228,7 @@ function assign_allocations!( # Save allocations to record push!(record.time, t) - push!(record.allocation_network_id, allocation_model.allocation_network_id) + push!(record.subnetwork_id, allocation_model.allocation_network_id) push!(record.user_node_id, Int(user_node_id)) push!(record.priority, allocation.priorities[priority_idx]) push!(record.demand, user.demand[user_idx]) @@ -503,7 +503,7 @@ function save_allocation_flows!( push!(record.edge_id, edge_metadata.id) push!(record.from_node_id, node_ids[i]) push!(record.to_node_id, node_ids[i + 1]) - push!(record.allocation_network_id, allocation_network_id) + push!(record.subnetwork_id, allocation_network_id) push!(record.priority, priority) push!(record.flow, flow) push!(record.collect_demands, collect_demands) @@ -518,7 +518,7 @@ function save_allocation_flows!( push!(record.edge_id, 0) push!(record.from_node_id, node_id) push!(record.to_node_id, node_id) - push!(record.allocation_network_id, allocation_network_id) + push!(record.subnetwork_id, allocation_network_id) push!(record.priority, priority) push!(record.flow, flow) push!(record.collect_demands, collect_demands) From dd273e60dbc9ad4fb98836471828a784fce375e2 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Fri, 16 Feb 2024 12:34:18 +0100 Subject: [PATCH 37/37] Fix test model --- core/src/read.jl | 4 ++++ python/ribasim_testmodels/ribasim_testmodels/allocation.py | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/core/src/read.jl b/core/src/read.jl index 990df32b7..85abc3651 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -205,6 +205,10 @@ function initialize_allocation!(p::Parameters, config::Config)::Nothing (; allocation_network_ids, allocation_models, main_network_connections) = allocation allocation_network_ids_ = sort(collect(keys(graph[].node_ids))) + if isempty(allocation_network_ids_) + return nothing + end + errors = non_positive_allocation_network_id(graph) if errors error("Allocation network initialization failed.") diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index ba5bc8107..285fa46e0 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -1662,7 +1662,7 @@ def allocation_target_model(): df=gpd.GeoDataFrame( data={ "type": node_type, - "allocation_network_id": 5 * [2], + "subnetwork_id": 5 * [2], }, index=pd.Index(np.arange(len(xy)) + 1, name="fid"), geometry=node_xy,