DiscreteControl

Set parameters of other nodes based on model state conditions (e.g. Basin level). The table below shows which parameters are controllable for a given node type.

using Ribasim
using DataFrames: DataFrame
using MarkdownTables
using PythonCall
# Render matplotlib figures off-screen; they are displayed by ending a cell with the figure
pyimport("matplotlib").use("Agg")
include("../../../utils/gen_docs_table.jl")

node_names = Symbol[]
controllable_parameters = String[]

for node_type in fieldtypes(Ribasim.ParametersIndependent)
    if node_type <: Ribasim.AbstractParameterNode
        node_name = nameof(node_type)
        controllable_fields = Ribasim.controllablefields(node_name)
        controllable_fields = sort!(string.(controllable_fields))
        if node_name == :TabulatedRatingCurve
            controllable_fields = map(
                s -> replace(s, "table" => "q(h) relationship (given by level, flow_rate)"),
                controllable_fields,
            )
        end
        if !isempty(controllable_fields)
            push!(node_names, Ribasim.snake_case(node_name))
            push!(controllable_parameters, join(controllable_fields, ", "))
        end
    end
end

df = DataFrame(:node => node_names, :controllable_parameters => controllable_parameters)

markdown_table(df)

node	controllable_parameters
linear_resistance	resistance
manning_resistance	manning_n
tabulated_rating_curve	q(h) relationship (given by level, flow_rate)
pump	flow_rate
outlet	flow_rate
pid_control	derivative, integral, proportional, target

1 Schematic

include("../../../utils/docs_plots.jl")

control_diagram = plot(;
    xlims = (0, 11), ylims = (-4, 7), size = (768, 730),
    aspect_ratio = :equal,
    axis = false, ticks = false, grid = false, legend = false,
    margin = 5Plots.mm,
)

control_states = ["off", "low", "high"]
active_state = "low"
state_flow_rates = [0.0, 0.5, 1.0]
selector_x, selector_y = 2.5, 4.5
pump_x = 8.35
basin_xs = [selector_x, pump_x]
basin_y = -2.2
radius = 1.6

for (i, state) in enumerate(control_states)
    active = state == active_state
    start_angle = -pi / 2 + (i - 1) * 2pi / length(control_states)
    end_angle = start_angle + 2pi / length(control_states)
    angles = range(start_angle, end_angle; length = 60)
    middle_angle = (start_angle + end_angle) / 2
    plot!(
        control_diagram,
        Shape(
            [selector_x; selector_x .+ radius .* cos.(angles); selector_x],
            [selector_y; selector_y .+ radius .* sin.(angles); selector_y],
        );
        fillcolor = active ? :lightgreen : :lightgrey,
        linecolor = :white, linewidth = 2,
    )
    annotate!(
        control_diagram,
        selector_x + 0.7radius * cos(middle_angle),
        selector_y + 0.7radius * sin(middle_angle),
        text(state, 11),
    )

    y = 6.1 - (i - 1) * 1.2
    plot!(
        control_diagram,
        Shape(
            [6.2, 10.5, 10.5, 6.2],
            [y - 0.45, y - 0.45, y + 0.45, y + 0.45],
        );
        fillcolor = active ? :lightgreen : :white,
        linecolor = active ? :darkgreen : :grey,
        linewidth = active ? 2 : 1,
    )
    annotate!(
        control_diagram, pump_x, y,
        text("$state: flow_rate = $(state_flow_rates[i])", 10),
    )
end

plot!(
    control_diagram, [selector_x, selector_x], [2.9, 1.8];
    color = :grey, linestyle = :dot,
)
plot!(
    control_diagram, [pump_x, pump_x], [3.25, 1.8];
    color = :grey, linestyle = :dot,
)
plot!(
    control_diagram, [selector_x, pump_x], [1.4, 1.4];
    color = "#808080", linewidth = 2,
)
for basin_x in basin_xs
    start_x = basin_x == selector_x ? basin_x : basin_x - 0.18
    start_y = basin_y + 0.15
    end_x = selector_x
    end_y = 1.4
    middle_x = (start_x + end_x) / 2
    middle_y = (start_y + end_y) / 2
    direction_x = end_x - start_x
    direction_y = end_y - start_y
    direction_length = hypot(direction_x, direction_y)
    unit_x = direction_x / direction_length
    unit_y = direction_y / direction_length
    perpendicular_x = -unit_y
    perpendicular_y = unit_x
    arrow_length = 0.16
    arrow_half_width = 0.07
    base_x = middle_x - unit_x * arrow_length / 2
    base_y = middle_y - unit_y * arrow_length / 2
    plot!(
        control_diagram,
        [start_x, end_x], [start_y, end_y];
        color = "#808080", linestyle = :dash, linewidth = 2,
    )
    plot!(
        control_diagram,
        Shape(
            [
                middle_x + unit_x * arrow_length / 2,
                base_x + perpendicular_x * arrow_half_width,
                base_x - perpendicular_x * arrow_half_width,
            ],
            [
                middle_y + unit_y * arrow_length / 2,
                base_y + perpendicular_y * arrow_half_width,
                base_y - perpendicular_y * arrow_half_width,
            ],
        );
        fillcolor = "#808080", linecolor = "#808080",
    )
end
scatter!(
    control_diagram, [(selector_x + pump_x) / 2], [1.4];
    markershape = :rtriangle, markersize = 6, color = "#808080",
)
scatter!(
    control_diagram, [selector_x], [1.4];
    markershape = :star5, markersize = 20, color = "#984EA3",
)
scatter!(
    control_diagram, [pump_x], [1.4];
    markershape = :hexagon, markersize = 18, color = :grey,
)
annotate!(control_diagram, 1.0, 1.4, text("DiscreteControl", 11))
annotate!(control_diagram, pump_x, 0.65, text("Pump", 11))
annotate!(
    control_diagram, (selector_x + pump_x) / 2, 1.8,
    text("control", 10),
)

for basin_x in basin_xs
    annotate!(
        control_diagram,
        basin_x == selector_x ? selector_x + 0.35 : 5.45,
        basin_x == selector_x ? -0.2 : -0.05,
        text("listen", 10),
    )
end

for basin_x in basin_xs
    plot!(
        control_diagram,
        Shape(
            [basin_x - 0.18, basin_x + 0.18, basin_x + 0.09, basin_x - 0.09],
            [basin_y + 0.15, basin_y + 0.15, basin_y - 0.15, basin_y - 0.15],
        );
        fillcolor = "#0072b2", linecolor = :black, linewidth = 1.5,
    )
    annotate!(
        control_diagram, basin_x, basin_y - 0.4, text("Basin", 10),
    )
end
control_diagram

2 Tables

2.1 Variable

The compound variable schema defines linear combinations of variables which can be used in conditions. This means that this schema defines new variables with the given compound_variable_id that look like \[ \text{weight}_1 * \text{variable}_1 + \text{weight}_2 * \text{variable}_2 + \ldots, \]

which can be for instance an average or a difference of variables. If a variable comes from a time-series, a look ahead \(\Delta t\) can be supplied.

table_from_schema(
    Ribasim.Schema.DiscreteControl.Variable;
    restrictions = Dict(
        :listen_node_id => "cannot be a Junction",
        :variable => "must be `level`, `storage`, `flow_rate`, `concentration.<substance>` or `concentration_external.<substance>`",
        :look_ahead => "only on transient boundary conditions, non-negative",
    ),
)

column	type	unit	optional	default	restriction
node_id	Int32	-	x	-	-
compound_variable_id	Int32	-	x	-	-
listen_node_id	Int32	-	x	-	cannot be a Junction
variable	String	-	x	-	must be level, storage, flow_rate, concentration.<substance> or concentration_external.<substance>
weight	Float64	-	✓	1.0	-
look_ahead	Float64	s	✓	0.0	only on transient boundary conditions, non-negative

These variables can be listened to:

• The storage of a Basin
• The level of a Basin
• The level of a LevelBoundary (supports look ahead)
• The flow rate through one of these node types: Pump, Outlet, TabulatedRatingCurve, LinearResistance, ManningResistance
• The inflow rate of a UserDemand
• The flow rate of a FlowBoundary (supports look ahead)
• The concentration of a substance in a Basin, with variable concentration.<substance>, e.g. concentration.Cl
• The external concentration of a substance in a Basin, with variable concentration_external.<substance>, as given in the Basin / concentration_external table

2.2 Condition

The condition schema defines conditions of the form the discrete_control node with this node_id listens to whether the variable given by the node_id and compound_variable_id is greater than threshold_high (at time \(t\)) when the condition was false previously, or threshold_low when true previously, adding hysteresis. Using threshold_low is optional and defaults to threshold_high. We recommend setting threshold_low below threshold_high to create a range that avoids frequent control state changes from small fluctuations and improves model performance. In equation form:

\[ \text{condition}=\begin{cases} \text{weight}_1 * \text{variable}_1 + \text{weight}_2 * \text{variable}_2 + \ldots > \text{greater\_than}(t), !\text{condition\_prev} \\ \text{weight}_1 * \text{variable}_1 + \text{weight}_2 * \text{variable}_2 + \ldots \geq \text{less\_than}(t), \text{condition\_prev} \end{cases} \]

Multiple conditions with different threshold_high and/or threshold_low values can be defined on the same compound_variable, each with their own unique condition_id.

At the start of the simulation, a condition is true if its variable exceeds the midpoint of threshold_low and threshold_high (which equals threshold_high without hysteresis).

table_from_schema(
    Ribasim.Schema.DiscreteControl.Condition;
    restrictions = Dict(:threshold_low => "≤ threshold_high"),
    defaults_override = Dict(:threshold_low => "threshold_high"),
)

column	type	unit	optional	default	restriction
node_id	Int32	-	x	-	-
compound_variable_id	Int32	-	x	-	-
condition_id	Int32	-	x	-	-
threshold_high	Float64	various	x	-	-
threshold_low	Float64	various	✓	threshold_high	≤ threshold_high
time	DateTime	-	✓	-	-

When no (or only a single) time is supplied for a given condition_id, the threshold_high value of this condition is unchanging over time. When multiple times are supplied, in between the given times the threshold_high and threshold_low values are block interpolated (forward fill), and outside they are constant given by the nearest time value. It is possible to configure a DiscreteControl node such that it very rapidly oscillates between control states, which isn’t realistic behavior. See the example model below.

from ribasim import Model, Node
from ribasim.nodes import basin, discrete_control, flow_boundary, pump
from shapely.geometry import Point
import matplotlib.pyplot as plt

model = Model(starttime="2020-01-01", endtime="2021-01-01", crs="EPSG:4326")

fb = model.flow_boundary.add(
    Node(1, Point(0, 0)), [flow_boundary.Static(flow_rate=[1.0])]
)

bsn = model.basin.add(
    Node(2, Point(1, 0)), [basin.Profile(area=[100.0, 100.0], level=[0.0, 1.0])]
)

pmp = model.pump.add(Node(3, Point(2, 0)), [pump.Static(flow_rate=[1.0])])

t = model.terminal.add(Node(4, Point(3, 0)))

dc = model.discrete_control.add(
    Node(5, Point(1.5, -1)),
    [
        discrete_control.Variable(
            compound_variable_id=[1], listen_node_id=[2], variable=["level"]
        )
    ],
)

model.link.add(fb, bsn)
model.link.add(bsn, pmp)
model.link.add(pmp, t)
model.link.add(dc, pmp)

ax = model.plot()
ax.figure.tight_layout()
ax.figure

Say there is constant inflow from the FlowBoundary node and there are no vertical forcings on the Basin. Say furthermore that the Pump is active only when the Basin level is above a reference value, with a flow greater than the boundary flow. In this case the Pump will rapidly turn on and off as the Basin level oscillates around the reference value.

To make you aware of these oscillations, Ribasim will give a warning if successive control state changes of a DiscreteControl node in simulated time are closer than min_discrete_control_interval, which is configurable and defaults to 1 second. The warning is given once per DiscreteControl node, for at most 10 nodes. Set min_discrete_control_interval to 0.0 to disable these warnings.

By default the conditions are evaluated after every timestep. To evaluate them only at fixed times, for instance daily, set the solver setting discrete_control_interval, see usage.

2.3 Logic

The logic schema defines which control states are triggered based on the truth of the conditions a DiscreteControl node listens to. DiscreteControl is applied in the Julia core as follows:

• During the simulation it is checked whether the truth of any of the conditions changes.
• When a condition changes, the corresponding DiscreteControl node ID is retrieved (node_id in the condition schema above).
• The truth value of all the conditions this DiscreteControl node listens to are retrieved, in the order of the condition IDs. This is then converted into a string of “T” for true and “F” for false. This string we call the truth state.
• The table below determines for the given DiscreteControl node ID and truth state what the corresponding control state is.
• For all the nodes this DiscreteControl node affects (as given by the “control” links in Links / static), their parameters are set to those parameters in NodeType / static corresponding to the determined control state.

table_from_schema(
    Ribasim.Schema.DiscreteControl.Logic;
    restrictions = Dict(
        :truth_state => "consists of the characters `T` (true), `F` (false) and `*` (any)",
    ),
)

column	type	optional	restriction
node_id	Int32	x	-
truth_state	String	x	consists of the characters T (true), F (false) and * (any)
control_state	String	x	-