From d1ed5500eb1c515eba23ccf0d4b6f27708ffa9e2 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 8 Aug 2023 09:41:04 +0200 Subject: [PATCH 01/12] solve.jl docstrings update --- core/src/solve.jl | 53 ++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 50 insertions(+), 3 deletions(-) diff --git a/core/src/solve.jl b/core/src/solve.jl index 4ebddea13..d56190c60 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -139,6 +139,12 @@ callback. Type parameter C indicates the content backing the StructVector, which can be a NamedTuple of Vectors or Arrow Primitives, and is added to avoid type instabilities. + +node_id: node ID of the TabulatedRatingCurve node +active: whether this node is active and thus contributes flows +tables: The current Q(h) relationships +time: The time table used for updating the tables +control_mapping: dictionary from (node_id, control_state) to Q(h) and/or active state """ struct TabulatedRatingCurve{C} <: AbstractParameterNode node_id::Vector{Int} @@ -153,6 +159,11 @@ Requirements: * from: must be (Basin,) node * to: must be (Basin,) node + +node_id: node ID of the LinearResistance node +active: whether this node is active and thus contributes flows +resistance: the resistance to flow; Q = Δh/resistance +control_mapping: dictionary from (node_id, control_state) to resistance and/or active state """ struct LinearResistance <: AbstractParameterNode node_id::Vector{Int} @@ -210,6 +221,10 @@ Requirements: * from: must be (TabulatedRatingCurve,) node * to: must be (Basin,) node * fraction must be positive. + +node_id: node ID of the TabulatedRatingCurve node +fraction: The fraction in [0,1] of flow the node lets through +control_mapping: dictionary from (node_id, control_state) to fraction """ struct FractionalFlow <: AbstractParameterNode node_id::Vector{Int} @@ -219,8 +234,8 @@ end """ node_id: node ID of the LevelBoundary node +active: whether this node is active level: the fixed level of this 'infinitely big basin' -The node_id are Indices to support fast lookup of level using ID. """ struct LevelBoundary <: AbstractParameterNode node_id::Vector{Int} @@ -230,8 +245,8 @@ end """ node_id: node ID of the FlowBoundary node +active: whether this node is active and thus contributes flow flow_rate: target flow rate -time: Data of time-dependent flow rates """ struct FlowBoundary <: AbstractParameterNode node_id::Vector{Int} @@ -241,8 +256,12 @@ end """ node_id: node ID of the Pump node +active: whether this node is active and thus contributes flow flow_rate: target flow rate +min_flow_rate: The minimal flow rate of the pump +max_flow_rate: The maximum flow rate of the pump control_mapping: dictionary from (node_id, control_state) to target flow rate +is_pid_controlled: whether the flow rate of this pump is governed by PID control """ struct Pump <: AbstractParameterNode node_id::Vector{Int} @@ -254,6 +273,25 @@ struct Pump <: AbstractParameterNode is_pid_controlled::BitVector end +""" +node_id: node ID of the Weir node +active: whether this node is active and thus contributes flow +flow_rate: target flow rate +min_flow_rate: The minimal flow rate of the weir +max_flow_rate: The maximum flow rate of the weir +control_mapping: dictionary from (node_id, control_state) to target flow rate +is_pid_controlled: whether the flow rate of this weir is governed by PID control +""" +struct Weir <: AbstractParameterNode + node_id::Vector{Int} + active::BitVector + flow_rate::Vector{Float64} + min_flow_rate::Vector{Float64} + max_flow_rate::Vector{Float64} + control_mapping::Dict{Tuple{Int, String}, NamedTuple} + is_pid_controlled:BitVector +end + """ node_id: node ID of the Terminal node """ @@ -262,7 +300,7 @@ struct Terminal <: AbstractParameterNode end """ -node_id: node ID of the Control node +node_id: node ID of the DiscreteControl node listen_feature_id: the ID of the node/edge being condition on variable: the name of the variable in the condition greater_than: The threshold value in the condition @@ -285,6 +323,15 @@ struct DiscreteControl <: AbstractParameterNode } end +""" +node_id: node ID of the PidControl node +active: whether this node is active and thus sets flow rates +listen_node_id: the id of the basin being controlled +proportional: the coefficient of the term proportional to the error +integral: the coefficient of the term proportional to the integral of the error +derivative: the coefficient of the term proportional to the derivative of the error +error: the current error; basin_target_level - current_level +""" struct PidControl <: AbstractParameterNode node_id::Vector{Int} active::BitVector From ce18802f0ae033e250eabecd179fbab34feaeaca Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 8 Aug 2023 16:58:08 +0200 Subject: [PATCH 02/12] Docs updates --- docs/contribute/addnode.qmd | 13 +++++++++++-- docs/core/equations.qmd | 22 +++++++++++----------- docs/core/usage.qmd | 6 ++++-- 3 files changed, 26 insertions(+), 15 deletions(-) diff --git a/docs/contribute/addnode.qmd b/docs/contribute/addnode.qmd index 59840749e..52cd86b29 100644 --- a/docs/contribute/addnode.qmd +++ b/docs/contribute/addnode.qmd @@ -29,6 +29,9 @@ There can be several schemas associated with a single node type. To define a sch ```julia @schema "ribasim.newnodetype.static" NewNodeTypeStatic +""" +node_id: node ID of the NewNodeType node +""" @version NewNodeTypeStaticV1 begin node_id::Int # Other fields @@ -44,9 +47,15 @@ Now we define the function that is called in the second bullet above, in `create ```julia function NewNodeType(db::DB, config::Config)::NewNodeType static = load_structvector(db, config, NewNodeTypeStaticV1) + defaults = (; active = true) + # Process potential control states in the static data + static_parsed = parse_static(static, db, "Weir", defaults) # Unpack the fields of static as inputs for the NewNodeType constructor - return NewNodeType(static.node_id, static.some_property) + return NewNodeType( + static_parsed.node_id, + static_parsed.some_property, + static_parsed.control_mapping) end ``` @@ -148,7 +157,7 @@ In `python/ribasim/ribasim/model.py`, add - `from ribasim.new_node_type import NewNodeType`; - new_node_type as a parameter and in the docstring of the `Model` class. -In `python/ribasim/ribasim/node.py` add a color and shape description in the `MARKERS` and `COLORS` dictionaries. +In `python/ribasim/ribasim/geometry/node.py` add a color and shape description in the `MARKERS` and `COLORS` dictionaries. # QGIS plugin diff --git a/docs/core/equations.qmd b/docs/core/equations.qmd index 8009dd789..7a8fce1e0 100644 --- a/docs/core/equations.qmd +++ b/docs/core/equations.qmd @@ -183,11 +183,11 @@ basins (external nodes) in a network. 2. `LevelBoundary` 3. `FlowBoundary` -### Pump +### Pump and weir -The behaviour of the pump is very straight forward if it is not PID controlled. Its flow is given by a fixed flow rate $q$, multiplied by a reduction factor: +The behaviour of pumps and weirs is very straight forward if these nodes are not PID controlled. Their flow is given by a fixed flow rate $q$, multiplied by a reduction factor: $$ -Q_\text{pump} = \phi(u)q +Q_\text{pump/weir} = \phi(u)q $$ Here $u$ is the storage of the upstream basin. The reduction factor @@ -415,10 +415,10 @@ $$ {#eq-PIDflow} for given constant parameters $K_p,K_i,K_d$. The pump or weir can have associated minimum and maximum flow rates $Q_\min, Q_\max$, and so $$ -Q_\text{pump/weir} = \text{clip}(\phi(u_\text{PID})Q_\text{PID}; Q_\min, Q_\max) +Q_\text{pump/weir} = \text{clip}(\phi(u_\text{us})Q_\text{PID}; Q_\min, Q_\max). $$ -where $u_\text{PID}$ is the storage of the controlled basin, $\phi$ is the reduction factor (see @eq-pumpreductionfactor) and +Here $u_\text{us}$ is the storage of the basin upstream of the pump or weir, $\phi$ is the reduction factor (see @eq-pumpreductionfactor) and $$ \text{clip}(Q; Q_\min, Q_\max) = @@ -450,16 +450,16 @@ $$ where $A(u_\text{PID})$ is the area of the controlled basin. The complexity arises from the fact that $Q_\text{PID}$ is a contribution to $\frac{\text{d}u_\text{PID}}{\text{d}t} = f_\text{PID}$, which makes @eq-PIDflow an implicit equation for $Q_\text{PID}$. We define $$ -f_\text{PID} = \hat{f}_\text{PID} - Q_\text{pump/weir}, +f_\text{PID} = \hat{f}_\text{PID} \pm Q_\text{pump/weir}, $$ -that is, $\hat{f}_\text{PID}$ is the right hand side of the ODE for the controlled basin storage state without the contribution of the PID controlled pump. The minus sign comes from the fact that the pump is puming away from the controlled basin. +that is, $\hat{f}_\text{PID}$ is the right hand side of the ODE for the controlled basin storage state without the contribution of the PID controlled pump. The plus sign holds for a weir and the minus sign for a pump, dictated by the way the pump and weir connectivity to the controlled basin is enforced. Using this, solving @eq-PIDflow for $Q_\text{PID}$ yields $$ -Q_\text{pump/weir} = \text{clip}\left(\phi(u_\text{PID})\frac{K_pe + K_iI + K_d \left(\frac{\text{d}\text{SP}}{\text{d}t}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right)}{1-\phi(u_\text{PID})\frac{K_d}{A(u_\text{PID})}};Q_\min,Q_\max\right), +Q_\text{pump/weir} = \text{clip}\left(\phi(u_\text{us})\frac{K_pe + K_iI + K_d \left(\frac{\text{d}\text{SP}}{\text{d}t}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right)}{1\pm\phi(u_\text{us})\frac{K_d}{A(u_\text{PID})}};Q_\min,Q_\max\right), $$ -where the clipping is again done last. Note that to compute this, $\hat{f}_\text{PID}$ has to be known first, meaning that the PID controlled pump flow rate has to be computed after all other contributions to the PID controlled basin's storage are known. +where the clipping is again done last. Note that to compute this, $\hat{f}_\text{PID}$ has to be known first, meaning that the PID controlled pump/weir flow rate has to be computed after all other contributions to the PID controlled basin's storage are known. ## Jacobian contributions @@ -475,7 +475,7 @@ Note that when the calculated pump flow lies outside $[Q_\min, Q_\max]$, the pum $$ \begin{align} E(u_\text{PID}) &= K_pe + K_iI + K_d \left(\dot{\text{SP}}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right) \\ - D(u_\text{PID}) &= 1-\phi(u_\text{PID})\frac{K_d}{A(u_\text{PID})} + D(u_\text{PID}) &= 1 \pm \phi(u_\text{PID})\frac{K_d}{A(u_\text{PID})} \end{align} $$ @@ -489,7 +489,7 @@ For the derivative with respect to $u_\text{PID}$ we need the derivative of the $$ \begin{align} E'(u_\text{PID}) &= -\frac{K_p}{A(u_\text{PID})} - K_d\frac{A(u_\text{PID})\frac{\partial\hat{f}_\text{PID}}{\partial u_\text{PID}} - \hat{f}_\text{PID}A'(u_\text{PID})}{A(u_\text{PID})^2}\\ - D'(u_\text{PID}) &= - K_d\left[\frac{\phi'(u_\text{PID})}{A(u_\text{PID})}- \phi(u_\text{PID})\frac{A'(u_\text{PID})}{A(u_\text{PID})^2}\right] + D'(u_\text{PID}) &= \pm K_d\left[\frac{\phi'(u_\text{PID})}{A(u_\text{PID})}- \phi(u_\text{PID})\frac{A'(u_\text{PID})}{A(u_\text{PID})^2}\right] \end{align} $$ For computational efficiency we exploit that $\phi'(u_\text{PID})$ is mostly $0$. Note that diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index c235d2093..93acf2bd7 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -276,9 +276,9 @@ level | Float64 | $m$ | - discharge | Float64 | $m^3 s^{-1}$ | non-negative -### Pump +### Pump and weir -Pump water from a source node to a destination node. +Pump/conduct water from a source node to a destination node. The set flow rate will be pumped unless the intake storage is less than $10~m^3$, in which case the flow rate will be linearly reduced to $0~m^3/s$. The intake must be either a Basin or LevelBoundary. @@ -459,6 +459,8 @@ control_state | String The PidControl node controls the level in a basin by continuously controlling the flow rate of a connected pump or weir. 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. For computational efficiency, if one of the terms should not be used, set its coefficient to `NaN` in stead of $0$. +The pump must be an outneighbor of the controlled basin, the weir must be an inneighbor. Thus for the same coefficient values in `proportional`, `integral` and `derivative`, the pump and weir have the oppositve effect on the controlled basin. + column | type | unit | restriction -------------- | -------- | -------- | ----------- node_id | Int | - | sorted From 399634cc0e60af86150f6ac236879ef38f095407 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 8 Aug 2023 17:04:55 +0200 Subject: [PATCH 03/12] Add weir to Ribasim python --- docs/schema/WeirStatic.schema.json | 57 +++++++++++++++++++ python/ribasim/ribasim/__init__.py | 2 + python/ribasim/ribasim/geometry/node.py | 2 + python/ribasim/ribasim/model.py | 4 ++ python/ribasim/ribasim/models.py | 50 +++++++++------- python/ribasim/ribasim/node_types/__init__.py | 2 + python/ribasim/ribasim/node_types/weir.py | 32 +++++++++++ 7 files changed, 129 insertions(+), 20 deletions(-) create mode 100644 docs/schema/WeirStatic.schema.json create mode 100644 python/ribasim/ribasim/node_types/weir.py diff --git a/docs/schema/WeirStatic.schema.json b/docs/schema/WeirStatic.schema.json new file mode 100644 index 000000000..d18f38d33 --- /dev/null +++ b/docs/schema/WeirStatic.schema.json @@ -0,0 +1,57 @@ +{ + "$schema": "https://json-schema.org/draft/2020-12/schema", + "properties": { + "max_flow_rate": { + "format": "default", + "description": "max_flow_rate", + "type": [ + "number" + ] + }, + "remarks": { + "format": "default", + "default": "", + "description": "a hack for pandera", + "type": "string" + }, + "active": { + "format": "default", + "description": "active", + "type": [ + "boolean" + ] + }, + "flow_rate": { + "format": "double", + "description": "flow_rate", + "type": "number" + }, + "node_id": { + "format": "default", + "description": "node_id", + "type": "integer" + }, + "control_state": { + "format": "default", + "description": "control_state", + "type": [ + "string" + ] + }, + "min_flow_rate": { + "format": "default", + "description": "min_flow_rate", + "type": [ + "number" + ] + } + }, + "required": [ + "node_id", + "flow_rate" + ], + "$id": "https://deltares.github.io/Ribasim/schema/WeirStatic.schema.json", + "title": "WeirStatic", + "description": "A WeirStatic object based on Ribasim.WeirStaticV1", + "type": "object" +} diff --git a/python/ribasim/ribasim/__init__.py b/python/ribasim/ribasim/__init__.py index e4128bb01..d3ddf13fd 100644 --- a/python/ribasim/ribasim/__init__.py +++ b/python/ribasim/ribasim/__init__.py @@ -16,6 +16,7 @@ from ribasim.node_types.pump import Pump from ribasim.node_types.tabulated_rating_curve import TabulatedRatingCurve from ribasim.node_types.terminal import Terminal +from ribasim.node_types.weir import Weir __all__ = [ "models", @@ -29,6 +30,7 @@ "Model", "Node", "Pump", + "Weir", "FlowBoundary", "Solver", "TabulatedRatingCurve", diff --git a/python/ribasim/ribasim/geometry/node.py b/python/ribasim/ribasim/geometry/node.py index 1410396c2..466dfe60c 100644 --- a/python/ribasim/ribasim/geometry/node.py +++ b/python/ribasim/ribasim/geometry/node.py @@ -105,6 +105,7 @@ def plot(self, ax=None, zorder=None) -> Any: "ManningResistance": "D", "TabulatedRatingCurve": "D", "Pump": "h", + "Weir": "h", "Terminal": "s", "FlowBoundary": "h", "DiscreteControl": "*", @@ -120,6 +121,7 @@ def plot(self, ax=None, zorder=None) -> Any: "ManningResistance": "r", "TabulatedRatingCurve": "g", "Pump": "0.5", # grayscale level + "Weir": "y", "Terminal": "m", "FlowBoundary": "m", "DiscreteControl": "k", diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py index 3b54ca055..756726f6c 100644 --- a/python/ribasim/ribasim/model.py +++ b/python/ribasim/ribasim/model.py @@ -28,6 +28,7 @@ from ribasim.node_types.pump import Pump from ribasim.node_types.tabulated_rating_curve import TabulatedRatingCurve from ribasim.node_types.terminal import Terminal +from ribasim.node_types.weir import Weir from ribasim.types import FilePath @@ -73,6 +74,8 @@ class Model(BaseModel): Tabulated rating curve describing flow based on the upstream water level. pump : Optional[Pump] Prescribed flow rate from one basin to the other. + weir : Optional[Weir] + Prescribed flow rate from one basin to the other. terminal : Optional[Terminal] Water sink without state or properties. discrete_control : Optional[DiscreteControl] @@ -98,6 +101,7 @@ class Model(BaseModel): manning_resistance: Optional[ManningResistance] tabulated_rating_curve: Optional[TabulatedRatingCurve] pump: Optional[Pump] + weir: Optional[Weir] terminal: Optional[Terminal] discrete_control: Optional[DiscreteControl] pid_control: Optional[PidControl] diff --git a/python/ribasim/ribasim/models.py b/python/ribasim/ribasim/models.py index 3aaa9843c..6ba1c50e8 100644 --- a/python/ribasim/ribasim/models.py +++ b/python/ribasim/ribasim/models.py @@ -1,6 +1,6 @@ # generated by datamodel-codegen: # filename: root.schema.json -# timestamp: 2023-07-20T13:32:54+00:00 +# timestamp: 2023-08-08T12:03:01+00:00 from __future__ import annotations @@ -42,6 +42,16 @@ class PumpStatic(BaseModel): min_flow_rate: Optional[float] = Field(None, description="min_flow_rate") +class WeirStatic(BaseModel): + max_flow_rate: Optional[float] = Field(None, description="max_flow_rate") + remarks: Optional[str] = Field("", description="a hack for pandera") + active: Optional[bool] = Field(None, description="active") + flow_rate: float = Field(..., description="flow_rate") + node_id: int = Field(..., description="node_id") + control_state: Optional[str] = Field(None, description="control_state") + min_flow_rate: Optional[float] = Field(None, description="min_flow_rate") + + class LevelBoundaryStatic(BaseModel): remarks: Optional[str] = Field("", description="a hack for pandera") active: Optional[bool] = Field(None, description="active") @@ -70,7 +80,6 @@ class BasinForcing(BaseModel): class FractionalFlowStatic(BaseModel): remarks: Optional[str] = Field("", description="a hack for pandera") - active: Optional[bool] = Field(None, description="active") node_id: int = Field(..., description="node_id") fraction: float = Field(..., description="fraction") control_state: Optional[str] = Field(None, description="control_state") @@ -166,22 +175,23 @@ class BasinStatic(BaseModel): class Root(BaseModel): - TerminalStatic: Optional[TerminalStatic] = None - TabulatedRatingCurveTime: Optional[TabulatedRatingCurveTime] = None - TabulatedRatingCurveStatic: Optional[TabulatedRatingCurveStatic] = None - PumpStatic: Optional[PumpStatic] = None - PidControlStatic: Optional[PidControlStatic] = None - Node: Optional[Node] = None - ManningResistanceStatic: Optional[ManningResistanceStatic] = None - LinearResistanceStatic: Optional[LinearResistanceStatic] = None - LevelBoundaryStatic: Optional[LevelBoundaryStatic] = None - FractionalFlowStatic: Optional[FractionalFlowStatic] = None - FlowBoundaryTime: Optional[FlowBoundaryTime] = None - FlowBoundaryStatic: Optional[FlowBoundaryStatic] = None - Edge: Optional[Edge] = None - DiscreteControlLogic: Optional[DiscreteControlLogic] = None - DiscreteControlCondition: Optional[DiscreteControlCondition] = None - BasinStatic: Optional[BasinStatic] = None - BasinState: Optional[BasinState] = None - BasinProfile: Optional[BasinProfile] = None BasinForcing: Optional[BasinForcing] = None + BasinProfile: Optional[BasinProfile] = None + BasinState: Optional[BasinState] = None + BasinStatic: Optional[BasinStatic] = None + DiscreteControlCondition: Optional[DiscreteControlCondition] = None + DiscreteControlLogic: Optional[DiscreteControlLogic] = None + Edge: Optional[Edge] = None + FlowBoundaryStatic: Optional[FlowBoundaryStatic] = None + FlowBoundaryTime: Optional[FlowBoundaryTime] = None + FractionalFlowStatic: Optional[FractionalFlowStatic] = None + LevelBoundaryStatic: Optional[LevelBoundaryStatic] = None + LinearResistanceStatic: Optional[LinearResistanceStatic] = None + ManningResistanceStatic: Optional[ManningResistanceStatic] = None + Node: Optional[Node] = None + PidControlStatic: Optional[PidControlStatic] = None + PumpStatic: Optional[PumpStatic] = None + TabulatedRatingCurveStatic: Optional[TabulatedRatingCurveStatic] = None + TabulatedRatingCurveTime: Optional[TabulatedRatingCurveTime] = None + TerminalStatic: Optional[TerminalStatic] = None + WeirStatic: Optional[WeirStatic] = None diff --git a/python/ribasim/ribasim/node_types/__init__.py b/python/ribasim/ribasim/node_types/__init__.py index cdf96049b..f5f91b6c2 100644 --- a/python/ribasim/ribasim/node_types/__init__.py +++ b/python/ribasim/ribasim/node_types/__init__.py @@ -9,6 +9,7 @@ from ribasim.node_types.pump import Pump from ribasim.node_types.tabulated_rating_curve import TabulatedRatingCurve from ribasim.node_types.terminal import Terminal +from ribasim.node_types.weir import Weir __all__ = [ "Basin", @@ -17,6 +18,7 @@ "LinearResistance", "ManningResistance", "Pump", + "Weir", "FlowBoundary", "TabulatedRatingCurve", "Terminal", diff --git a/python/ribasim/ribasim/node_types/weir.py b/python/ribasim/ribasim/node_types/weir.py new file mode 100644 index 000000000..276029553 --- /dev/null +++ b/python/ribasim/ribasim/node_types/weir.py @@ -0,0 +1,32 @@ +import pandera as pa +from pandera.engines.pandas_engine import PydanticModel +from pandera.typing import DataFrame + +from ribasim import models +from ribasim.input_base import TableModel + +__all__ = ("Weir",) + + +class StaticSchema(pa.SchemaModel): + class Config: + """Config with dataframe-level data type.""" + + dtype = PydanticModel(models.WeirStatic) + + +class Weir(TableModel): + """ + Conducts water from a source node to a destination node. + The set flow rate will be used unless the intake storage is less than 10m3, + in which case the flow rate will be linearly reduced to 0 m3/s. + Negative flow rates are not supported. + Note that the intake must always be a Basin. + + Parameters + ---------- + static : pd.DataFrame + Table with constant flow rates. + """ + + static: DataFrame[StaticSchema] From c816949989c3a559a2fd9ec88075f9c0a75eaf6a Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 8 Aug 2023 17:05:30 +0200 Subject: [PATCH 04/12] QGIS update --- qgis/core/nodes.py | 36 ++++++++++++++++++++++++++++++++++-- 1 file changed, 34 insertions(+), 2 deletions(-) diff --git a/qgis/core/nodes.py b/qgis/core/nodes.py index d483024f8..5b4676d88 100644 --- a/qgis/core/nodes.py +++ b/qgis/core/nodes.py @@ -177,6 +177,7 @@ def renderer(self) -> QgsCategorizedSymbolRenderer: "LevelBoundary": (QColor("green"), "LevelBoundary", shape.Circle), "FlowBoundary": (QColor("purple"), "FlowBoundary", shape.Hexagon), "Pump": (QColor("gray"), "Pump", shape.Hexagon), + "Weir": (QColor("yellow"), "Weir", shape.Hexagon), "ManningResistance": (QColor("red"), "ManningResistance", shape.Diamond), "Terminal": (QColor("purple"), "Terminal", shape.Square), "DiscreteControl": (QColor("black"), "DiscreteControl", shape.Star), @@ -325,8 +326,10 @@ class TabulatedRatingCurveStatic(Input): geometry_type = "No Geometry" attributes = [ QgsField("node_id", QVariant.Int), + QgsField("active", QVariant.Bool), QgsField("level", QVariant.Double), QgsField("discharge", QVariant.Double), + QgsField("control_state", QVariant.String), ] @@ -356,6 +359,7 @@ class LinearResistanceStatic(Input): geometry_type = "No Geometry" attributes = [ QgsField("node_id", QVariant.Int), + QgsField("active", QVariant.Bool), QgsField("resistance", QVariant.Double), QgsField("control_state", QVariant.String), ] @@ -366,6 +370,7 @@ class ManningResistanceStatic(Input): geometry_type = "No Geometry" attributes = [ QgsField("node_id", QVariant.Int), + QgsField("active", QVariant.Bool), QgsField("length", QVariant.Double), QgsField("manning_n", QVariant.Double), QgsField("profile_width", QVariant.Double), @@ -379,8 +384,8 @@ class LevelBoundaryStatic(Input): geometry_type = "No Geometry" attributes = [ QgsField("node_id", QVariant.Int), + QgsField("active", QVariant.Bool), QgsField("level", QVariant.Double), - QgsField("control_state", QVariant.String), ] @@ -389,7 +394,23 @@ class PumpStatic(Input): geometry_type = "No Geometry" attributes = [ QgsField("node_id", QVariant.Int), + QgsField("active", QVariant.Bool), + QgsField("flow_rate", QVariant.Double), + QgsField("min_flow_rate", QVariant.Double), + QgsField("max_flow_rate", QVariant.Double), + QgsField("control_state", QVariant.String), + ] + + +class WeirStatic(Input): + input_type = "Weir / static" + geometry_type = "No Geometry" + attributes = [ + QgsField("node_id", QVariant.Int), + QgsField("active", QVariant.Bool), QgsField("flow_rate", QVariant.Double), + QgsField("min_flow_rate", QVariant.Double), + QgsField("max_flow_rate", QVariant.Double), QgsField("control_state", QVariant.String), ] @@ -404,9 +425,19 @@ class FlowBoundaryStatic(Input): input_type = "FlowBoundary / static" geometry_type = "No Geometry" attributes = [ + QgsField("node_id", QVariant.Int), + QgsField("active", QVariant.Bool), + QgsField("flow_rate", QVariant.Double), + ] + + +class FlowBoundaryTime(Input): + input_type = "FlowBoundary / time" + geometry_type = "No Geometry" + attributes = [ + QgsField("time", QVariant.DateTime), QgsField("node_id", QVariant.Int), QgsField("flow_rate", QVariant.Double), - QgsField("control_state", QVariant.String), ] @@ -436,6 +467,7 @@ class PidControlStatic(Input): geometry_type = "LineString" attributes = [ QgsField("node_id", QVariant.Int), + QgsField("active", QVariant.Bool), QgsField("listen_node_id", QVariant.Int), QgsField("proportional", QVariant.Double), QgsField("integral", QVariant.Double), From ef5a490540af13c3737ee9a6980f7c283f2d1a36 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 8 Aug 2023 17:05:48 +0200 Subject: [PATCH 05/12] PID control testmodel update --- .../ribasim_testmodels/pid_control.py | 55 ++++++++++++------- 1 file changed, 35 insertions(+), 20 deletions(-) diff --git a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py index ad118c8ec..26beb493e 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py +++ b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py @@ -11,9 +11,11 @@ def pid_control_model(): [ (0.0, 0.0), # 1: FlowBoundary (1.0, 0.0), # 2: Basin - (2.0, 0.0), # 3: Pump + (2.0, 0.5), # 3: Pump (3.0, 0.0), # 4: LevelBoundary (1.5, 1.0), # 5: PidControl + (2.0, -0.5), # 6: Weir + (1.5, -1.0), # 7: PidControl ] ) @@ -25,6 +27,8 @@ def pid_control_model(): "Pump", "LevelBoundary", "PidControl", + "Weir", + "PidControl", ] # Make sure the feature id starts at 1: explicitly give an index. @@ -38,8 +42,8 @@ def pid_control_model(): ) # Setup the edges: - from_id = np.array([1, 2, 3, 5], dtype=np.int64) - to_id = np.array([2, 3, 4, 3], dtype=np.int64) + from_id = np.array([1, 2, 3, 4, 6, 5, 7], dtype=np.int64) + to_id = np.array([2, 3, 4, 6, 2, 3, 6], dtype=np.int64) lines = ribasim.utils.geometry_from_connectivity(node, from_id, to_id) edge = ribasim.Edge( @@ -47,7 +51,7 @@ def pid_control_model(): data={ "from_node_id": from_id, "to_node_id": to_id, - "edge_type": 3 * ["flow"] + ["control"], + "edge_type": 5 * ["flow"] + 2 * ["control"], }, geometry=lines, crs="EPSG:28992", @@ -55,14 +59,9 @@ def pid_control_model(): ) # Setup the basins: - R = 10.0 - n = 10 - - level = np.linspace(0, R, n) - area = np.pi * level * (2 * R - level) - area[0] = 0.01 - - profile = pd.DataFrame(data={"node_id": n * [2], "level": level, "area": area}) + profile = pd.DataFrame( + data={"node_id": [2, 2], "level": [0.0, 1.0], "area": [1000.0, 1000.0]} + ) # Convert steady forcing to m/s # 2 mm/d precipitation, 1 mm/d evaporation @@ -78,11 +77,18 @@ def pid_control_model(): "infiltration": [0.0], "precipitation": [precipitation], "urban_runoff": [0.0], - "target_level": [R / 2], + "target_level": [5.0], + } + ) + + state = pd.DataFrame( + data={ + "node_id": [2], + "storage": [500.0], } ) - basin = ribasim.Basin(profile=profile, static=static) + basin = ribasim.Basin(profile=profile, static=static, state=state) # Setup pump: pump = ribasim.Pump( @@ -90,7 +96,16 @@ def pid_control_model(): data={ "node_id": [3], "flow_rate": [0.0], # Will be overwritten by PID controller - "min_flow_rate": [0.0], + } + ) + ) + + # Setup weir: + weir = ribasim.Weir( + static=pd.DataFrame( + data={ + "node_id": [6], + "flow_rate": [0.0], # Will be overwritten by PID controller } ) ) @@ -114,11 +129,10 @@ def pid_control_model(): pid_control = ribasim.PidControl( static=pd.DataFrame( data={ - "node_id": [5], - "listen_node_id": [2], - "proportional": [-1e-3], - "derivative": [None], - "integral": [-1e-7], + "node_id": [5, 7], + "listen_node_id": [2, 2], + "proportional": [-1e-3, 1e-3], + "integral": [-1e-7, 1e-7], } ) ) @@ -132,6 +146,7 @@ def pid_control_model(): flow_boundary=flow_boundary, level_boundary=level_boundary, pump=pump, + weir=weir, pid_control=pid_control, starttime="2020-01-01 00:00:00", endtime="2020-07-01 00:00:00", From 47466886b78423b3034bb65374aeba40c4d2e128 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 8 Aug 2023 17:06:42 +0200 Subject: [PATCH 06/12] Small improvements --- core/test/validation.jl | 4 +++- docs/schema/FractionalFlowStatic.schema.json | 7 ------- docs/schema/root.schema.json | 3 +++ python/ribasim/ribasim/node_types/pump.py | 2 +- 4 files changed, 7 insertions(+), 9 deletions(-) diff --git a/core/test/validation.jl b/core/test/validation.jl index 3481253ed..e068f958f 100644 --- a/core/test/validation.jl +++ b/core/test/validation.jl @@ -88,6 +88,7 @@ end @testset "PidControl connectivity validation" begin pid_control_node_id = [1, 6] pid_control_listen_node_id = [3, 5] + pump_node_id = [2, 4] graph_flow = DiGraph(7) graph_control = DiGraph(7) @@ -108,6 +109,7 @@ end graph_flow, graph_control, basin_node_id, + pump_node_id, ) end @@ -116,7 +118,7 @@ end @test logger.logs[1].message == "Listen node #3 of PidControl node #1 is not a Basin" @test logger.logs[2].level == Error @test logger.logs[2].message == - "Listen node #5 of PidControl node #6 is not upstream of controlled node #2" + "Listen node #5 of PidControl node #6 is not upstream of controlled pump #2" end # This test model is not written on Ubuntu CI, see #479 diff --git a/docs/schema/FractionalFlowStatic.schema.json b/docs/schema/FractionalFlowStatic.schema.json index d1e80f695..f98567832 100644 --- a/docs/schema/FractionalFlowStatic.schema.json +++ b/docs/schema/FractionalFlowStatic.schema.json @@ -7,13 +7,6 @@ "description": "a hack for pandera", "type": "string" }, - "active": { - "format": "default", - "description": "active", - "type": [ - "boolean" - ] - }, "node_id": { "format": "default", "description": "node_id", diff --git a/docs/schema/root.schema.json b/docs/schema/root.schema.json index b27ea9070..19db7c6c0 100644 --- a/docs/schema/root.schema.json +++ b/docs/schema/root.schema.json @@ -13,6 +13,9 @@ "PumpStatic": { "$ref": "PumpStatic.schema.json" }, + "WeirStatic": { + "$ref": "WeirStatic.schema.json" + }, "LevelBoundaryStatic": { "$ref": "LevelBoundaryStatic.schema.json" }, diff --git a/python/ribasim/ribasim/node_types/pump.py b/python/ribasim/ribasim/node_types/pump.py index cd169970d..a6ad5ff56 100644 --- a/python/ribasim/ribasim/node_types/pump.py +++ b/python/ribasim/ribasim/node_types/pump.py @@ -20,7 +20,7 @@ class Pump(TableModel): Pump water from a source node to a destination node. The set flow rate will be pumped unless the intake storage is less than 10m3, in which case the flow rate will be linearly reduced to 0 m3/s. - A negative flow rate means pumping against the edge direction. + Negative flow rates are not supported. Note that the intake must always be a Basin. Parameters From e76900ba506784b67ea32b11fb5e47f2e5623e64 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 8 Aug 2023 17:07:19 +0200 Subject: [PATCH 07/12] Weir implementation without jac --- core/src/bmi.jl | 17 +++++++--- core/src/create.jl | 29 +++++++++++++++-- core/src/jac.jl | 28 ++++++++-------- core/src/solve.jl | 74 ++++++++++++++++++++++++++++-------------- core/src/utils.jl | 52 +++++++++++++++++++++-------- core/src/validation.jl | 42 +++++++++++++++++++----- core/test/control.jl | 28 +++++++++++++--- 7 files changed, 196 insertions(+), 74 deletions(-) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index dda0c0137..5aced8d95 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -21,13 +21,14 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model error("Invalid number of connections for certain node types.") end - (; pid_control, connectivity, basin, pump, fractional_flow) = parameters + (; pid_control, connectivity, basin, pump, weir, fractional_flow) = parameters if !valid_pid_connectivity( pid_control.node_id, pid_control.listen_node_id, connectivity.graph_flow, connectivity.graph_control, basin.node_id, + pump.node_id, ) error("Invalid PidControl connectivity.") end @@ -41,9 +42,14 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model end for id in pid_control.node_id - id_pump = only(outneighbors(connectivity.graph_control, id)) - pump_idx = findsorted(pump.node_id, id_pump) - pump.is_pid_controlled[pump_idx] = true + id_controlled = only(outneighbors(connectivity.graph_control, id)) + pump_idx = findsorted(pump.node_id, id_controlled) + if isnothing(pump_idx) + weir_idx = findsorted(weir.node_id, id_controlled) + weir.is_pid_controlled[weir_idx] = true + else + pump.is_pid_controlled[pump_idx] = true + end end # tstops for transient flow_boundary @@ -74,7 +80,8 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model timespan = (zero(t_end), t_end) jac_prototype = get_jac_prototype(parameters) - RHS = ODEFunction(water_balance!; jac_prototype, jac = water_balance_jac!) + # TODO: Undo this and update jac.jl + RHS = ODEFunction(water_balance!; jac_prototype)#, jac = water_balance_jac!) @timeit_debug to "Setup ODEProblem" begin prob = ODEProblem(RHS, u0, timespan, parameters) diff --git a/core/src/create.jl b/core/src/create.jl index 0d6b0be2d..38d63e6be 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -4,11 +4,19 @@ function parse_static( nodetype::String, defaults::NamedTuple, )::NamedTuple + # E.g. `PumpStatic` static_type = eltype(static) columnnames_static = collect(fieldnames(static_type)) mask = [symb ∉ [:node_id, :control_state] for symb in columnnames_static] + + # The names of the variables that can define a control state columnnames_variables = columnnames_static[mask] + + # The types of the variables that can define a control state columntypes_variables = collect(fieldtypes(static_type))[mask] + + # A vector of vectors, for each variable the (initial) values for all nodes + # of the current type vals = [] node_ids = get_ids(db, nodetype) @@ -273,8 +281,6 @@ function Pump(db::DB, config::Config)::Pump static = load_structvector(db, config, PumpStaticV1) defaults = (; min_flow_rate = 0.0, max_flow_rate = NaN, active = true) static_parsed = parse_static(static, db, "Pump", defaults) - - # TODO: use this in formulate_jac! for pump is_pid_controlled = falses(length(static_parsed.node_id)) return Pump( @@ -288,6 +294,23 @@ function Pump(db::DB, config::Config)::Pump ) end +function Weir(db::DB, config::Config)::Weir + static = load_structvector(db, config, WeirStaticV1) + defaults = (; min_flow_rate = 0.0, max_flow_rate = NaN, active = true) + static_parsed = parse_static(static, db, "Weir", defaults) + is_pid_controlled = falses(length(static_parsed.node_id)) + + return Weir( + static_parsed.node_id, + static_parsed.active, + static_parsed.flow_rate, + static_parsed.min_flow_rate, + static_parsed.max_flow_rate, + static_parsed.control_mapping, + is_pid_controlled, + ) +end + function Terminal(db::DB, config::Config)::Terminal static = load_structvector(db, config, TerminalStaticV1) return Terminal(static.node_id) @@ -407,6 +430,7 @@ function Parameters(db::DB, config::Config)::Parameters level_boundary = LevelBoundary(db, config) flow_boundary = FlowBoundary(db, config) pump = Pump(db, config) + weir = Weir(db, config) terminal = Terminal(db, config) discrete_control = DiscreteControl(db, config) pid_control = PidControl(db, config) @@ -424,6 +448,7 @@ function Parameters(db::DB, config::Config)::Parameters level_boundary, flow_boundary, pump, + weir, terminal, discrete_control, pid_control, diff --git a/core/src/jac.jl b/core/src/jac.jl index 971836961..d0181d6f4 100644 --- a/core/src/jac.jl +++ b/core/src/jac.jl @@ -212,17 +212,17 @@ function formulate_jac!( end """ -The contributions of Pump nodes to the Jacobian. +The contributions of Pump and Weir nodes to the Jacobian. """ function formulate_jac!( - pump::Pump, + node::Union{Pump, Weir}, J::SparseMatrixCSC{Float64, Int64}, u::ComponentVector{Float64}, p::Parameters, t::Float64, )::Nothing - (; basin, fractional_flow, connectivity, pid_control) = p - (; active, node_id, flow_rate, is_pid_controlled) = pump + (; basin, fractional_flow, connectivity) = p + (; active, node_id, flow_rate, is_pid_controlled) = node (; graph_flow) = connectivity @@ -357,7 +357,7 @@ function formulate_jac!( end """ -The contributions of TabulatedRatingCurve nodes to the Jacobian. +The contributions of PidControl nodes to the Jacobian. """ function formulate_jac!( pid_control::PidControl, @@ -366,8 +366,7 @@ function formulate_jac!( p::Parameters, t::Float64, )::Nothing - # TODO: update after pid_control equation test merge - (; basin, connectivity, pump) = p + (; basin, connectivity, pump, weir) = p (; node_id, active, listen_node_id, proportional, integral, derivative, error) = pid_control (; min_flow_rate, max_flow_rate) = pump @@ -391,18 +390,17 @@ function formulate_jac!( continue end - id_pump = only(outneighbors(graph_control, id)) + id_controlled = only(outneighbors(graph_control, id)) listened_node_id = listen_node_id[i] _, listened_node_idx = id_index(basin.node_id, listened_node_id) listen_area = basin.current_area[listened_node_idx] storage_listened_basin = u.storage[listened_node_idx] - area = basin.current_area[listened_node_idx] reduction_factor = min(storage_listened_basin, 10.0) / 10.0 K_d = derivative[i] if !isnan(K_d) # dlevel/dstorage = 1/area - D = 1.0 - K_d * reduction_factor / area + D = 1.0 - K_d * reduction_factor / listen_area else D = 1.0 end @@ -422,7 +420,7 @@ function formulate_jac!( if !isnan(K_d) dtarget_level = 0.0 du_listened_basin_old = du.storage[listened_node_idx] - E += K_d * (dtarget_level - du_listened_basin_old / area) + E += K_d * (dtarget_level - du_listened_basin_old / listen_area) end # Clip values outside pump flow rate bounds @@ -461,23 +459,23 @@ function formulate_jac!( dD = reduction_factor * darea if reduction_factor_regime - dD -= dreduction_factor / area + dD -= dreduction_factor / listen_area end dD *= K_d dE = -K_d * ( - area * J[listened_node_idx, listened_node_idx] - + listen_area * J[listened_node_idx, listened_node_idx] - du.storage[listened_node_idx] * darea - ) / (area^2) + ) / (listen_area^2) else dD = 0.0 dE = 0.0 end if !isnan(K_p) - dE -= K_p / area + dE -= K_p / listen_area end if reduction_factor_regime diff --git a/core/src/solve.jl b/core/src/solve.jl index d56190c60..6e407c487 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -289,7 +289,7 @@ struct Weir <: AbstractParameterNode min_flow_rate::Vector{Float64} max_flow_rate::Vector{Float64} control_mapping::Dict{Tuple{Int, String}, NamedTuple} - is_pid_controlled:BitVector + is_pid_controlled::BitVector end """ @@ -354,6 +354,7 @@ struct Parameters level_boundary::LevelBoundary flow_boundary::FlowBoundary pump::Pump + weir::Weir terminal::Terminal discrete_control::DiscreteControl pid_control::PidControl @@ -519,7 +520,7 @@ function continuous_control!( )::Nothing # TODO: Also support being able to control weir # TODO: also support time varying target levels - (; connectivity, pump, basin, fractional_flow) = p + (; connectivity, pump, weir, basin, fractional_flow) = p (; min_flow_rate, max_flow_rate) = pump (; graph_control, graph_flow, flow) = connectivity (; node_id, active, proportional, integral, derivative, listen_node_id, error) = @@ -528,16 +529,21 @@ function continuous_control!( get_error!(pid_control, p) for (i, id) in enumerate(node_id) - controlled_node_id = only(outneighbors(graph_control, id)) - # TODO: support the use of id_index - controlled_node_idx = searchsortedfirst(pump.node_id, controlled_node_id) - if !active[i] du.integral[i] = 0.0 u.integral[i] = 0.0 return end + controlled_node_id = only(outneighbors(graph_control, id)) + controls_pump = (controlled_node_id in pump.node_id) + + controlled_node_idx = if controls_pump + findsorted(pump.node_id, controlled_node_id) + else + findsorted(weir.node_id, controlled_node_id) + end + du.integral[i] = error[i] listened_node_id = listen_node_id[i] @@ -576,14 +582,22 @@ function continuous_control!( end # Clip values outside pump flow rate bounds - flow_rate = clamp(flow_rate, min_flow_rate[i], max_flow_rate[i]) + flow_rate = clamp( + flow_rate, + min_flow_rate[controlled_node_idx], + max_flow_rate[controlled_node_idx], + ) # Below du.storage is updated. This is normally only done # in formulate!(du, connectivity, basin), but in this function # flows are set so du has to be updated too. - - pump.flow_rate[controlled_node_idx] = flow_rate - du.storage[listened_node_idx] -= flow_rate + if controls_pump + pump.flow_rate[controlled_node_idx] = flow_rate + du.storage[listened_node_idx] -= flow_rate + else + weir.flow_rate[controlled_node_idx] = flow_rate + du.storage[listened_node_idx] += flow_rate + end # Set flow for connected edges src_id = only(inneighbors(graph_flow, controlled_node_id)) @@ -597,17 +611,20 @@ function continuous_control!( du.storage[dst_idx] += flow_rate end - for id in outneighbors(graph_flow, controlled_node_id) - if id in fractional_flow.node_id - after_ff_id = only(outneighbours(graph_flow, id)) - ff_idx = findsorted(fractional_flow, id) - flow_rate_fraction = fractional_flow.fraction[ff_idx] * flow_rate - flow[id, after_ff_id] = flow_rate_fraction + # When the controlled pump flows out into fractional flow nodes + if controls_pump + for id in outneighbors(graph_flow, controlled_node_id) + if id in fractional_flow.node_id + after_ff_id = only(outneighbours(graph_flow, id)) + ff_idx = findsorted(fractional_flow, id) + flow_rate_fraction = fractional_flow.fraction[ff_idx] * flow_rate + flow[id, after_ff_id] = flow_rate_fraction - has_index, basin_idx = id_index(basin.node_id, after_ff_id) + has_index, basin_idx = id_index(basin.node_id, after_ff_id) - if has_index - du.storage[basin_idx] += flow_rate_fraction + if has_index + du.storage[basin_idx] += flow_rate_fraction + end end end end @@ -786,13 +803,21 @@ function formulate!(flow_boundary::FlowBoundary, p::Parameters, t::Float64)::Not end end -function formulate!(pump::Pump, p::Parameters, storage::AbstractVector{Float64})::Nothing - (; connectivity, basin, level_boundary) = p +function formulate!( + node::Union{Pump, Weir}, + p::Parameters, + storage::AbstractVector{Float64}, +)::Nothing + (; connectivity, basin) = p (; graph_flow, flow) = connectivity - (; node_id, active, flow_rate, is_pid_controlled) = pump + (; node_id, active, flow_rate, is_pid_controlled) = node for (id, isactive, rate, pid_controlled) in zip(node_id, active, flow_rate, is_pid_controlled) - @assert rate >= 0 "Pump flow rate must be positive, found $rate for Pump #$id" + if rate < 0 + error( + "$(typeof(node)) flow rate must be non-negative, found $rate for Pump #$id.", + ) + end src_id = only(inneighbors(graph_flow, id)) dst_id = only(outneighbors(graph_flow, id)) @@ -816,7 +841,6 @@ function formulate!(pump::Pump, p::Parameters, storage::AbstractVector{Float64}) q = reduction_factor * rate else # Pumping from level boundary - @assert src_id in level_boundary.node_id "Pump intake is neither basin nor level_boundary" q = rate end @@ -858,6 +882,7 @@ function formulate_flows!( fractional_flow, flow_boundary, pump, + weir, ) = p formulate!(linear_resistance, p) @@ -866,6 +891,7 @@ function formulate_flows!( formulate!(flow_boundary, p, t) formulate!(fractional_flow, p) formulate!(pump, p, storage) + formulate!(weir, p, storage) return nothing end diff --git a/core/src/utils.jl b/core/src/utils.jl index 2d502193c..98a3790dd 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -528,14 +528,18 @@ Upstream basins always depend on themselves. function update_jac_prototype!( jac_prototype::SparseMatrixCSC{Float64, Int64}, p::Parameters, - node::Union{Pump, TabulatedRatingCurve}, + node::Union{Pump, Weir, TabulatedRatingCurve}, )::Nothing (; basin, fractional_flow, connectivity) = p (; graph_flow) = connectivity - for id in node.node_id + for (i, id) in enumerate(node.node_id) id_in = only(inneighbors(graph_flow, id)) + if hasfield(typeof(node), :is_pid_controlled) && node.is_pid_controlled[i] + continue + end + # For inneighbors only directly connected basins give a contribution has_index_in, idx_in = id_index(basin.node_id, id_in) @@ -576,14 +580,17 @@ function update_jac_prototype!( p::Parameters, node::PidControl, )::Nothing - (; basin, connectivity) = p + (; basin, connectivity, pump) = p (; graph_control, graph_flow) = connectivity n_basins = length(basin.node_id) - for (pid_idx, (listen_node_id, id)) in enumerate(zip(node.listen_node_id, node.node_id)) - id_pump = only(outneighbors(graph_control, id)) - id_pump_out = only(inneighbors(graph_flow, id_pump)) + for i in eachindex(node.node_id) + listen_node_id = node.listen_node_id[i] + id = node.node_id[i] + + # ID of controlled pump/weir + id_controlled = only(outneighbors(graph_control, id)) _, listen_idx = id_index(basin.node_id, listen_node_id) @@ -591,19 +598,36 @@ function update_jac_prototype!( jac_prototype[listen_idx, listen_idx] = 1.0 # PID control integral state - pid_state_idx = n_basins + pid_idx + pid_state_idx = n_basins + i jac_prototype[listen_idx, pid_state_idx] = 1.0 jac_prototype[pid_state_idx, listen_idx] = 1.0 - # The basin downstream of the pump - has_index, idx_out_out = id_index(basin.node_id, id_pump_out) + if id_controlled in pump.node_id + id_pump_out = only(inneighbors(graph_flow, id_controlled)) + + # The basin downstream of the pump + has_index, idx_out_out = id_index(basin.node_id, id_pump_out) + + if has_index + # The basin downstream of the pump depends on PID control integral state + jac_prototype[pid_state_idx, idx_out_out] = 1.0 + + # The basin downstream of the pump also depends on the controlled basin + jac_prototype[listen_idx, idx_out_out] = 1.0 + end + else + id_weir_in = only(outneighbors(graph_flow, id_controlled)) - if has_index - # The basin downstream of the pump PID control integral state - jac_prototype[pid_state_idx, idx_out_out] = 1.0 + # The basin upstream of the weir + has_index, idx_out_in = id_index(basin.node_id, id_weir_in) - # The basin downstream of the pump also depends on the controlled basin - jac_prototype[listen_idx, idx_out_out] = 1.0 + if has_index + # The basin upstream of the weir depends on the PID control integral state + jac_prototype[pid_state_idx, idx_out_in] = 1.0 + + # The basin upstream of the weir also depends on the controlled basin + jac_prototype[listen_idx, idx_out_in] = 1.0 + end end end return nothing diff --git a/core/src/validation.jl b/core/src/validation.jl index c2497e47a..4707cda30 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -4,7 +4,6 @@ @schema "ribasim.edge" Edge @schema "ribasim.discretecontrol.condition" DiscreteControlCondition @schema "ribasim.discretecontrol.logic" DiscreteControlLogic -@schema "ribasim.pump.static" PumpStatic @schema "ribasim.basin.static" BasinStatic @schema "ribasim.basin.forcing" BasinForcing @schema "ribasim.basin.profile" BasinProfile @@ -17,8 +16,10 @@ @schema "ribasim.linearresistance.static" LinearResistanceStatic @schema "ribasim.manningresistance.static" ManningResistanceStatic @schema "ribasim.pidcontrol.static" PidControlStatic +@schema "ribasim.pump.static" PumpStatic @schema "ribasim.tabulatedratingcurve.static" TabulatedRatingCurveStatic @schema "ribasim.tabulatedratingcurve.time" TabulatedRatingCurveTime +@schema "ribasim.weir.static" WeirStatic const delimiter = " / " tablename(sv::Type{SchemaVersion{T, N}}) where {T, N} = join(nodetype(sv), delimiter) @@ -42,23 +43,26 @@ end # Allowed types for downstream (to_node_id) nodes given the type of the upstream (from_node_id) node neighbortypes(nodetype::Symbol) = neighbortypes(Val(nodetype)) neighbortypes(::Val{:Pump}) = Set((:Basin, :FractionalFlow, :Terminal, :LevelBoundary)) +neighbortypes(::Val{:Weir}) = Set((:Basin, :FractionalFlow, :Terminal, :LevelBoundary)) neighbortypes(::Val{:Basin}) = - Set((:LinearResistance, :TabulatedRatingCurve, :ManningResistance, :Pump)) + Set((:LinearResistance, :TabulatedRatingCurve, :ManningResistance, :Pump, :Weir)) neighbortypes(::Val{:Terminal}) = Set{Symbol}() # only endnode neighbortypes(::Val{:FractionalFlow}) = Set((:Basin, :Terminal, :LevelBoundary)) neighbortypes(::Val{:FlowBoundary}) = Set((:Basin, :FractionalFlow, :Terminal, :LevelBoundary)) -neighbortypes(::Val{:LevelBoundary}) = Set((:LinearResistance, :ManningResistance, :Pump)) +neighbortypes(::Val{:LevelBoundary}) = + Set((:LinearResistance, :ManningResistance, :Pump, :Weir)) neighbortypes(::Val{:LinearResistance}) = Set((:Basin, :LevelBoundary)) neighbortypes(::Val{:ManningResistance}) = Set((:Basin, :LevelBoundary)) neighbortypes(::Val{:DiscreteControl}) = Set(( :Pump, + :Weir, :TabulatedRatingCurve, :LinearResistance, :ManningResistance, :FractionalFlow, )) -neighbortypes(::Val{:PidControl}) = Set((:Pump,)) +neighbortypes(::Val{:PidControl}) = Set((:Pump, :Weir)) neighbortypes(::Val{:TabulatedRatingCurve}) = Set((:Basin, :FractionalFlow, :Terminal, :LevelBoundary)) neighbortypes(::Any) = Set{Symbol}() @@ -82,6 +86,7 @@ n_neighbor_bounds(::Val{:LevelBoundary}) = n_neighbor_bounds(::Val{:FlowBoundary}) = n_neighbor_bounds(0, 0, 1, typemax(Int)) neighbourtypes(::Any) = n_neighbor_bounds(0, 0, 0, 0) n_neighbor_bounds(::Val{:Pump}) = n_neighbor_bounds(1, 1, 1, typemax(Int)) +n_neighbor_bounds(::Val{:Weir}) = n_neighbor_bounds(1, 1, 1, typemax(Int)) n_neighbor_bounds(::Val{:Terminal}) = n_neighbor_bounds(1, typemax(Int), 0, 0) n_neighbor_bounds(::Val{:PidControl}) = n_neighbor_bounds(0, 0, 1, 1) n_neighbor_bounds(::Val{:DiscreteControl}) = n_neighbor_bounds(0, 0, 1, typemax(Int)) @@ -108,6 +113,15 @@ end control_state::Union{Missing, String} end +@version WeirStaticV1 begin + node_id::Int + active::Union{Missing, Bool} + flow_rate::Float64 + min_flow_rate::Union{Missing, Float64} + max_flow_rate::Union{Missing, Float64} + control_state::Union{Missing, String} +end + @version BasinStaticV1 begin node_id::Int drainage::Float64 @@ -351,21 +365,31 @@ function valid_pid_connectivity( graph_flow::DiGraph{Int}, graph_control::DiGraph{Int}, basin_node_id::Indices{Int}, + pump_node_id::Vector{Int}, )::Bool errors = false for (id, listen_id) in zip(pid_control_node_id, pid_control_listen_node_id) - pump_id = only(outneighbors(graph_control, id)) has_index, _ = id_index(basin_node_id, listen_id) if !has_index @error "Listen node #$listen_id of PidControl node #$id is not a Basin" errors = true end - pump_intake_id = only(inneighbors(graph_flow, pump_id)) - if pump_intake_id != listen_id - @error "Listen node #$listen_id of PidControl node #$id is not upstream of controlled node #$pump_id" - errors = true + controlled_id = only(outneighbors(graph_control, id)) + + if controlled_id in pump_node_id + pump_intake_id = only(inneighbors(graph_flow, controlled_id)) + if pump_intake_id != listen_id + @error "Listen node #$listen_id of PidControl node #$id is not upstream of controlled pump #$controlled_id" + errors = true + end + else + weir_outflow_id = only(outneighbors(graph_flow, controlled_id)) + if weir_outflow_id != listen_id + @error "Listen node #$listen_id of PidControl node #$id is not downstream of controlled weir #$controlled_id" + errors = true + end end end diff --git a/core/test/control.jl b/core/test/control.jl index b1cc56d65..2f50dca7b 100644 --- a/core/test/control.jl +++ b/core/test/control.jl @@ -57,12 +57,30 @@ end toml_path = normpath(@__DIR__, "../../data/pid_control/pid_control.toml") @test ispath(toml_path) model = Ribasim.run(toml_path) - basin = model.integrator.p.basin + p = model.integrator.p + (; basin, pid_control, flow_boundary) = p + + storage = Ribasim.get_storages_and_levels(model).storage[1, :] + timesteps = Ribasim.timesteps(model) + + K_p = pid_control.proportional[2] + K_i = pid_control.integral[2] + A = basin.area[1][1] + target_level = basin.target_level[1] + initial_storage = storage[1] + flow_rate = flow_boundary.flow_rate[1].u[1] + du0 = flow_rate + K_p * (target_level - initial_storage / A) + target_storage = A * target_level + Δstorage = initial_storage - target_storage + alpha = -K_p / (2 * A) + omega = sqrt(4 * K_i / A - (K_i / A)^2) / 2 + phi = atan(du0 / Δstorage - alpha) / omega + a = abs(Δstorage / cos(phi)) + bound = @. a * exp(alpha * timesteps) + # TODO: Make smaller after jac.jl update + eps = 3.0 - timesteps = Ribasim.timesteps(model) / (60 * 60 * 24) - level = Ribasim.get_storages_and_levels(model).level[1, :] - bound = 5 .* exp.(-0.03 .* timesteps) - @test all(abs.(level .- basin.target_level[1]) .< bound) + @test all((storage .- target_storage) .< bound .+ eps) end @testset "TabulatedRatingCurve control" begin From f094798fffcf0ae5f1961c5254538795de413876 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 9 Aug 2023 13:17:27 +0200 Subject: [PATCH 08/12] Update Jacobian --- core/src/bmi.jl | 3 +- core/src/jac.jl | 117 +++++++++++++----- core/src/solve.jl | 31 +++-- core/test/control.jl | 2 +- docs/core/equations.qmd | 60 +++++++-- .../ribasim_testmodels/pid_control.py | 10 +- 6 files changed, 161 insertions(+), 62 deletions(-) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 5aced8d95..cd3ac9e7b 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -80,8 +80,7 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model timespan = (zero(t_end), t_end) jac_prototype = get_jac_prototype(parameters) - # TODO: Undo this and update jac.jl - RHS = ODEFunction(water_balance!; jac_prototype)#, jac = water_balance_jac!) + RHS = ODEFunction(water_balance!; jac_prototype, jac = water_balance_jac!) @timeit_debug to "Setup ODEProblem" begin prob = ODEProblem(RHS, u0, timespan, parameters) diff --git a/core/src/jac.jl b/core/src/jac.jl index d0181d6f4..6e987689d 100644 --- a/core/src/jac.jl +++ b/core/src/jac.jl @@ -390,17 +390,40 @@ function formulate_jac!( continue end - id_controlled = only(outneighbors(graph_control, id)) + controlled_node_id = only(outneighbors(graph_control, id)) + controls_pump = (controlled_node_id in pump.node_id) + listened_node_id = listen_node_id[i] _, listened_node_idx = id_index(basin.node_id, listened_node_id) listen_area = basin.current_area[listened_node_idx] - storage_listened_basin = u.storage[listened_node_idx] - reduction_factor = min(storage_listened_basin, 10.0) / 10.0 + + if controls_pump + controlled_node_idx = findsorted(pump.node_id, controlled_node_id) + + listened_basin_storage = u.storage[listened_node_idx] + reduction_factor = min(listened_basin_storage, 10.0) / 10.0 + else + controlled_node_idx = findsorted(weir.node_id, controlled_node_id) + + # Upstream node of weir does not have to be a basin + upstream_node_id = only(inneighbors(graph_flow, controlled_node_id)) + has_upstream_index, upstream_basin_idx = + id_index(basin.node_id, upstream_node_id) + if has_upstream_index + upstream_basin_storage = u.storage[upstream_basin_idx] + reduction_factor = min(upstream_basin_storage, 10.0) / 10.0 + else + reduction_factor = 1.0 + end + end K_d = derivative[i] if !isnan(K_d) - # dlevel/dstorage = 1/area - D = 1.0 - K_d * reduction_factor / listen_area + if controls_pump + D = 1.0 - K_d * reduction_factor / listen_area + else + D = 1.0 + K_d * reduction_factor / listen_area + end else D = 1.0 end @@ -427,15 +450,15 @@ function formulate_jac!( flow_rate = reduction_factor * E / D was_clipped = false - if flow_rate < min_flow_rate[i] + if flow_rate < min_flow_rate[controlled_node_idx] was_clipped = true - flow_rate = min_flow_rate[i] + flow_rate = min_flow_rate[controlled_node_idx] end - if !isnan(max_flow_rate[i]) - if flow_rate > max_flow_rate[i] + if !isnan(max_flow_rate[controlled_node_idx]) + if flow_rate > max_flow_rate[controlled_node_idx] was_clipped = true - flow_rate = max_flow_rate[i] + flow_rate = max_flow_rate[controlled_node_idx] end end @@ -451,51 +474,87 @@ function formulate_jac!( end # Only in this case the reduction factor has a non-zero derivative - reduction_factor_regime = (storage_listened_basin < 10.0) + reduction_factor_regime = if controls_pump + listened_basin_storage < 10.0 + else + if has_upstream_index + upstream_basin_storage < 10.0 + else + false + end + end # Computing D and E derivatives if !isnan(K_d) darea = basin.current_darea[listened_node_idx] - dD = reduction_factor * darea + + dD_du_listen = reduction_factor * darea / (listen_area^2) if reduction_factor_regime - dD -= dreduction_factor / listen_area + if controls_pump + dD_du_listen -= 0.1 / darea + else + dD_du_upstream = -0.1 * K_d / area + end + else + dD_du_upstream = 0.0 end - dD *= K_d + dD_du_listen *= K_d - dE = + dE_du_listen = -K_d * ( listen_area * J[listened_node_idx, listened_node_idx] - du.storage[listened_node_idx] * darea ) / (listen_area^2) else - dD = 0.0 - dE = 0.0 + dD_du_listen = 0.0 + dD_du_upstream = 0.0 + dE_du_listen = 0.0 end if !isnan(K_p) - dE -= K_p / listen_area + dE_du_listen -= K_p / listen_area + end + + dQ_du_listen = reduction_factor * (D * dE_du_listen - E * dD_du_listen) / (D^2) + + if controls_pump && reduction_factor_regime + dQ_du_listen += 0.1 * E / D end - if reduction_factor_regime - dreduction_factor = 0.1 + if controls_pump + J[listened_node_idx, listened_node_idx] -= dQ_du_listen - dq = dreduction_factor * E / D + downstream_node_id = only(outneighbors(graph_flow, controlled_node_id)) + has_downstream_index, downstream_node_idx = + id_index(basin.node_id, downstream_node_id) + + if has_downstream_index + J[listened_node_idx, downstream_node_idx] += dQ_du_listen + end else - dq = 0.0 + J[listened_node_idx, listened_node_idx] += dQ_du_listen + + if has_upstream_index + J[listened_node_idx, upstream_basin_idx] -= dQ_du_listen + end end - dq += reduction_factor * (D * dE - E * dD) / (D^2) + if !controls_pump + if !isnan(K_d) && has_upstream_index + dE_du_upstream = -K_d * J[upstream_basin_idx, listened_node_idx] / area - J[listened_node_idx, listened_node_idx] -= dq + dQ_du_upstream = + reduction_factor * (D * dE_du_upstream - E * dD_du_upstream) / (D^2) - id_out = only(outneighbors(graph_flow, id_pump)) - has_index, idx_out_out = id_index(basin.node_id, id_out) + if reduction_factor_regime + dQ_du_upstream += 0.1 * E / D + end - if has_index - J[pid_state_idx, idx_out_out] += K_i * reduction_factor / D - J[listened_node_idx, idx_out_out] += dq + J[upstream_basin_idx, listened_node_idx] += dQ_du_upstream + J[upstream_basin_idx, upstream_basin_idx] -= dQ_du_upstream + end end end return nothing diff --git a/core/src/solve.jl b/core/src/solve.jl index 6e407c487..8d91321d0 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -535,22 +535,33 @@ function continuous_control!( return end + du.integral[i] = error[i] + + listened_node_id = listen_node_id[i] + _, listened_node_idx = id_index(basin.node_id, listened_node_id) + controlled_node_id = only(outneighbors(graph_control, id)) controls_pump = (controlled_node_id in pump.node_id) - controlled_node_idx = if controls_pump - findsorted(pump.node_id, controlled_node_id) + if controls_pump + controlled_node_idx = findsorted(pump.node_id, controlled_node_id) + + listened_basin_storage = u.storage[listened_node_idx] + reduction_factor = min(listened_basin_storage, 10.0) / 10.0 else - findsorted(weir.node_id, controlled_node_id) + controlled_node_idx = findsorted(weir.node_id, controlled_node_id) + + # Upstream node of weir does not have to be a basin + upstream_node_id = only(inneighbors(graph_flow, controlled_node_id)) + has_index, upstream_basin_idx = id_index(basin.node_id, upstream_node_id) + if has_index + upstream_basin_storage = u.storage[upstream_basin_idx] + reduction_factor = min(upstream_basin_storage, 10.0) / 10.0 + else + reduction_factor = 1.0 + end end - du.integral[i] = error[i] - - listened_node_id = listen_node_id[i] - _, listened_node_idx = id_index(basin.node_id, listened_node_id) - storage_listened_basin = u.storage[listened_node_idx] - reduction_factor = min(storage_listened_basin, 10.0) / 10.0 - flow_rate = 0.0 K_d = derivative[i] diff --git a/core/test/control.jl b/core/test/control.jl index 2f50dca7b..0fece78d4 100644 --- a/core/test/control.jl +++ b/core/test/control.jl @@ -76,8 +76,8 @@ end omega = sqrt(4 * K_i / A - (K_i / A)^2) / 2 phi = atan(du0 / Δstorage - alpha) / omega a = abs(Δstorage / cos(phi)) + # This bound is the exact envelope of the analytical solution bound = @. a * exp(alpha * timesteps) - # TODO: Make smaller after jac.jl update eps = 3.0 @test all((storage .- target_storage) .< bound .+ eps) diff --git a/docs/core/equations.qmd b/docs/core/equations.qmd index 7a8fce1e0..13945236f 100644 --- a/docs/core/equations.qmd +++ b/docs/core/equations.qmd @@ -441,13 +441,17 @@ $$ \frac{\text{d}I}{\text{d}t} = e(t). $$ +::: {.callout-note} +In the case of the controlled weir, the upstream node can also be a level boundary. In this case we define $\phi = 1$. +::: + ## The derivative term When $K_d \ne 0$ this adds a level of complexity. We can see this by looking at the error derivative more closely: $$ \frac{\text{d}e}{\text{d}t} = \frac{\text{d}\text{SP}}{\text{d}t} - \frac{1}{A(u_\text{PID})}\frac{\text{d}u_\text{PID}}{\text{d}t}, $$ -where $A(u_\text{PID})$ is the area of the controlled basin. The complexity arises from the fact that $Q_\text{PID}$ is a contribution to $\frac{\text{d}u_\text{PID}}{\text{d}t} = f_\text{PID}$, which makes @eq-PIDflow an implicit equation for $Q_\text{PID}$. We define +where $A(u_\text{PID})$ is the area of the controlled basin as a function of the storage of the controlled basin $u_\text{PID}$. The complexity arises from the fact that $Q_\text{PID}$ is a contribution to $\frac{\text{d}u_\text{PID}}{\text{d}t} = f_\text{PID}$, which makes @eq-PIDflow an implicit equation for $Q_\text{PID}$. We define $$ f_\text{PID} = \hat{f}_\text{PID} \pm Q_\text{pump/weir}, @@ -474,34 +478,66 @@ Note that when the calculated pump flow lies outside $[Q_\min, Q_\max]$, the pum $$ \begin{align} - E(u_\text{PID}) &= K_pe + K_iI + K_d \left(\dot{\text{SP}}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right) \\ - D(u_\text{PID}) &= 1 \pm \phi(u_\text{PID})\frac{K_d}{A(u_\text{PID})} + E &= K_pe + K_iI + K_d \left(\dot{\text{SP}}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right) \\ + D &= 1 \pm \phi(u_\text{us})\frac{K_d}{A(u_\text{PID})} \end{align} $$ Then $$ -\frac{\partial Q_\text{pump/weir}}{\partial I} = K_i\frac{\phi(u_\text{PID})}{D(u_\text{PID})}. +\frac{\partial Q_\text{pump/weir}}{\partial I} = K_i\frac{\phi(u_\text{us})}{D}. $$ -For the derivative with respect to $u_\text{PID}$ we need the derivative of the enumerator and denominator: +For the derivative of $Q_\text{pump/weir}$ with respect to storages we need the derivatives of the enumerator and denominator. + +For the denominator we distinguish 2 possibilities: + +- If the controlled node is a pump, then $u_\text{us} \equiv u_\text{PID}$ and so $D$ is only a function of $u_\text{PID}$: +$$ +\frac{\text{d}D}{\text{d}u_\text{PID}} = - K_d\left[\frac{\phi'(u_\text{PID})}{A(u_\text{PID})}- \phi(u_\text{PID})\frac{A'(u_\text{PID})}{A(u_\text{PID})^2}\right] +$$ +- If the controlled node is a weir, then $D$ is a function of two distinct storages (assuming the upstream node is indeed a basin): $$ \begin{align} - E'(u_\text{PID}) &= -\frac{K_p}{A(u_\text{PID})} - K_d\frac{A(u_\text{PID})\frac{\partial\hat{f}_\text{PID}}{\partial u_\text{PID}} - \hat{f}_\text{PID}A'(u_\text{PID})}{A(u_\text{PID})^2}\\ - D'(u_\text{PID}) &= \pm K_d\left[\frac{\phi'(u_\text{PID})}{A(u_\text{PID})}- \phi(u_\text{PID})\frac{A'(u_\text{PID})}{A(u_\text{PID})^2}\right] + \frac{\partial D}{\partial u_\text{PID}} &= \phi(u_\text{us})K_d\frac{A'(u_\text{PID})}{A(u_\text{PID})^2}, \\ + \frac{\partial D}{\partial u_\text{us}} &= -\phi'(u_\text{us})K_d\frac{1}{A(u_\text{PID})}. \end{align} $$ -For computational efficiency we exploit that $\phi'(u_\text{PID})$ is mostly $0$. Note that +For the enumerator there distinguish 2 cases: + +- The partial derivative with respect to $u_\text{PID}$: +$$ +\frac{\partial E}{\partial u_\text{PID}} = -\frac{K_p}{A(u_\text{PID})} - K_d\frac{A(u_\text{PID})\frac{\partial\hat{f}_\text{PID}}{\partial u_\text{PID}} - \hat{f}_\text{PID}A'(u_\text{PID})}{A(u_\text{PID})^2}. $$ - \frac{\partial\hat{f}_\text{PID}}{\partial u_\text{PID}} = \hat{J}[\text{PID},\text{PID}], +- The partial derivative with respect to a different storage $u_i$: $$ -\emph{i.e.} the Jacobian term of the dependence of $u_\text{PID}$ on itself without the contribution of the PID controlled pump/weir. Note thus that also the Jacobian contribution of the PID controlled pump can only be computed after the other contributions to the Jacobian of $u_\text{PID}$ to itself are known. +\frac{\partial E}{\partial u_i} = - \frac{K_d}{A(u_\text{PID})}\frac{\partial \hat{f}_\text{PID}}{\partial u_i}. +$$ + +For computational efficiency we exploit that $\phi'$ is mostly $0$. Note that -The partial derivative with respect to $u_\text{PID}$ is then given by $$ -\frac{\partial Q_\text{pump/weir}}{\partial u_\text{PID}} = \phi'(u_\text{PID})\frac{E(u_\text{PID})}{D(u_\text{PID})} + \phi(u_\text{PID})\frac{D(u_\text{PID})E'(u_\text{PID}) - E(u_\text{PID})D'(u_\text{PID})}{D(u_\text{PID})^2}. +\frac{\partial\hat{f}_\text{PID}}{\partial u_i} = \hat{J}[i,\text{PID}], +$$ +\emph{i.e.} the Jacobian term of the dependence of $u_\text{PID}$ on $u_i$ without the contribution of the PID controlled pump/weir. Note thus that also the Jacobian contribution of the PID controlled pump can only be computed after the other contributions to the Jacobian are known. + +The partial derivatives of $Q_\text{pump}$ are then given by +$$ +\begin{align} + \frac{\partial Q_\text{pump}}{\partial u_\text{PID}} &= \phi'(u_\text{PID})\frac{E}{D} + \frac{\phi(u_\text{PID})}{D^2}\left[D\frac{\partial E}{\partial u_\text{PID}} - E\frac{\text{d}D}{\text{d}u_\text{PID}}\right], \\ + \frac{\partial Q_\text{pump}}{\partial u_i} &= \frac{\phi(u_\text{PID})}{D}\frac{\partial E}{\partial u_i} \qquad (i \ne \text{PID}). +\end{align} +$$ + +The partial derivatives of $Q_\text{weir}$ are given by +$$ +\begin{align} + \frac{\partial Q_\text{weir}}{\partial u_\text{PID}} &= \frac{\phi(u_\text{us})}{D^2}\left[D\frac{\partial E}{\partial u_\text{PID}} - E\frac{\partial D}{\partial u_\text{PID}}\right], \\ + \frac{\partial Q_\text{weir}}{\partial u_\text{us}} &= \phi'(u_\text{us})\frac{E}{D} + \frac{\phi(u_\text{us})}{D^2}\left[D\frac{\partial E}{\partial u_\text{us}} - E\frac{\partial D}{\partial u_\text{us}}\right], \\ + \frac{\partial Q_\text{weir}}{\partial u_i} &= \frac{\phi(u_\text{us})}{D}\frac{\partial E}{\partial u_i} \qquad (i \ne \text{PID, us}). +\end{align} $$ ## The sign of the parameters diff --git a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py index 26beb493e..b7844858d 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py +++ b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py @@ -63,19 +63,13 @@ def pid_control_model(): data={"node_id": [2, 2], "level": [0.0, 1.0], "area": [1000.0, 1000.0]} ) - # Convert steady forcing to m/s - # 2 mm/d precipitation, 1 mm/d evaporation - seconds_in_day = 24 * 3600 - precipitation = 0.002 / seconds_in_day - evaporation = 0.001 / seconds_in_day - static = pd.DataFrame( data={ "node_id": [2], "drainage": [0.0], - "potential_evaporation": [evaporation], + "potential_evaporation": [0.0], "infiltration": [0.0], - "precipitation": [precipitation], + "precipitation": [0.0], "urban_runoff": [0.0], "target_level": [5.0], } From 6c141c57e5aa7c083a4173fb6b231fec41b0038e Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 9 Aug 2023 13:20:13 +0200 Subject: [PATCH 09/12] Fix tests --- core/test/utils.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/core/test/utils.jl b/core/test/utils.jl index 4492a5200..fbc7800fc 100644 --- a/core/test/utils.jl +++ b/core/test/utils.jl @@ -126,9 +126,9 @@ end p = Ribasim.Parameters(db, cfg) jac_prototype = Ribasim.get_jac_prototype(p) - @test jac_prototype.m == 2 - @test jac_prototype.n == 2 - @test jac_prototype.colptr == [1, 3, 4] - @test jac_prototype.rowval == [1, 2, 1] - @test jac_prototype.nzval == ones(3) + @test jac_prototype.m == 3 + @test jac_prototype.n == 3 + @test jac_prototype.colptr == [1, 4, 5, 6] + @test jac_prototype.rowval == [1, 2, 3, 1, 1] + @test jac_prototype.nzval == ones(5) end From 1073c34d87783ba477e608b87e4feddeac030fe3 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 10 Aug 2023 16:29:50 +0200 Subject: [PATCH 10/12] Bring back valid_discrete_control --- core/src/bmi.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index cd3ac9e7b..44bbd3d26 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -21,6 +21,10 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model error("Invalid number of connections for certain node types.") end + if !valid_discrete_control(parameters) + error("Invalid discrete control state definition(s).") + end + (; pid_control, connectivity, basin, pump, weir, fractional_flow) = parameters if !valid_pid_connectivity( pid_control.node_id, From 062701639bf79daa2b156328ec7245c035e8582f Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 10 Aug 2023 17:20:34 +0200 Subject: [PATCH 11/12] Refactor flow_rate validation --- core/src/solve.jl | 54 +++++++++++++++++++++++++++++++++++++----- core/src/validation.jl | 40 +++++++++++++++++++++++++++++++ 2 files changed, 88 insertions(+), 6 deletions(-) diff --git a/core/src/solve.jl b/core/src/solve.jl index 8d91321d0..51398e2b9 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -271,6 +271,30 @@ struct Pump <: AbstractParameterNode max_flow_rate::Vector{Float64} control_mapping::Dict{Tuple{Int, String}, NamedTuple} is_pid_controlled::BitVector + + function Pump( + node_id, + active, + flow_rate, + min_flow_rate, + max_flow_rate, + control_mapping, + is_pid_controlled, + ) + if valid_flow_rates(node_id, flow_rate, control_mapping, :Pump) + return new( + node_id, + active, + flow_rate, + min_flow_rate, + max_flow_rate, + control_mapping, + is_pid_controlled, + ) + else + error("Invalid Pump flow rate(s).") + end + end end """ @@ -290,6 +314,30 @@ struct Weir <: AbstractParameterNode max_flow_rate::Vector{Float64} control_mapping::Dict{Tuple{Int, String}, NamedTuple} is_pid_controlled::BitVector + + function Weir( + node_id, + active, + flow_rate, + min_flow_rate, + max_flow_rate, + control_mapping, + is_pid_controlled, + ) + if valid_flow_rates(node_id, flow_rate, control_mapping, :Weir) + return new( + node_id, + active, + flow_rate, + min_flow_rate, + max_flow_rate, + control_mapping, + is_pid_controlled, + ) + else + error("Invalid Weir flow rate(s).") + end + end end """ @@ -824,12 +872,6 @@ function formulate!( (; node_id, active, flow_rate, is_pid_controlled) = node for (id, isactive, rate, pid_controlled) in zip(node_id, active, flow_rate, is_pid_controlled) - if rate < 0 - error( - "$(typeof(node)) flow rate must be non-negative, found $rate for Pump #$id.", - ) - end - src_id = only(inneighbors(graph_flow, id)) dst_id = only(outneighbors(graph_flow, id)) diff --git a/core/src/validation.jl b/core/src/validation.jl index 4707cda30..edf82d1d3 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -359,6 +359,46 @@ function valid_profiles( return errors end +""" +Test whether static or discrete controlled flow rates are indeed non-negative. +""" +function valid_flow_rates( + node_id::Vector{Int}, + flow_rate::Vector{Float64}, + control_mapping::Dict{Tuple{Int, String}, NamedTuple}, + node_type::Symbol, +)::Bool + errors = false + + # Collect ids of discrete controlled nodes so that they do not give another error + # if their initial value is also invalid. + ids_controlled = Int[] + + for (key, control_values) in pairs(control_mapping) + id_controlled = key[1] + push!(ids_controlled, id_controlled) + flow_rate_ = get(control_values, :flow_rate, 1) + + if flow_rate_ < 0.0 + errors = true + control_state = key[2] + @error "Flow rates must be non-negative, found $flow_rate_ for control state '$control_state' of $node_type #$id_controlled." + end + end + + for (id, flow_rate_) in zip(node_id, flow_rate) + if id in ids_controlled + continue + end + if flow_rate_ < 0.0 + errors = true + @error "Flow rates must be non-negative, found $flow_rate_ for static $node_type #$id." + end + end + + return !errors +end + function valid_pid_connectivity( pid_control_node_id::Vector{Int}, pid_control_listen_node_id::Vector{Int}, From a2319c870d6ac89a3bddfc5bdee1b6ea69fc04b0 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 10 Aug 2023 18:23:59 +0200 Subject: [PATCH 12/12] Add flow_rate sign validation test --- core/src/validation.jl | 4 ++-- core/test/validation.jl | 41 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+), 2 deletions(-) diff --git a/core/src/validation.jl b/core/src/validation.jl index edf82d1d3..2e16205c4 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -382,7 +382,7 @@ function valid_flow_rates( if flow_rate_ < 0.0 errors = true control_state = key[2] - @error "Flow rates must be non-negative, found $flow_rate_ for control state '$control_state' of $node_type #$id_controlled." + @error "$node_type flow rates must be non-negative, found $flow_rate_ for control state '$control_state' of #$id_controlled." end end @@ -392,7 +392,7 @@ function valid_flow_rates( end if flow_rate_ < 0.0 errors = true - @error "Flow rates must be non-negative, found $flow_rate_ for static $node_type #$id." + @error "$node_type flow rates must be non-negative, found $flow_rate_ for static #$id." end end diff --git a/core/test/validation.jl b/core/test/validation.jl index 801404691..e2531bd28 100644 --- a/core/test/validation.jl +++ b/core/test/validation.jl @@ -178,3 +178,44 @@ end @test logger.logs[1].message == "These control states from DiscreteControl node #4 are not defined for controlled Ribasim.Pump #2: [\"foo\"]." end + +@testset "Pump/weir flow rate sign validation" begin + logger = TestLogger() + + with_logger(logger) do + @test_throws "Invalid Weir flow rate(s)." Ribasim.Weir( + [1], + [true], + [-1.0], + [NaN], + [NaN], + Dict{Tuple{Int, String}, NamedTuple}(), + [false], + ) + end + + @assert length(logger.logs) == 1 + @test logger.logs[1].level == Error + @test logger.logs[1].message == + "Weir flow rates must be non-negative, found -1.0 for static #1." + + logger = TestLogger() + + with_logger(logger) do + @test_throws "Invalid Pump flow rate(s)." Ribasim.Pump( + [1], + [true], + [-1.0], + [NaN], + [NaN], + Dict{Tuple{Int, String}, NamedTuple}((1, "foo") => (; flow_rate = -1.0)), + [false], + ) + end + + # Only the invalid control state flow_rate yields an error + @assert length(logger.logs) == 1 + @test logger.logs[1].level == Error + @test logger.logs[1].message == + "Pump flow rates must be non-negative, found -1.0 for control state 'foo' of #1." +end