diff --git a/core/src/callback.jl b/core/src/callback.jl index 242c1eafa..22c5fbb49 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -9,10 +9,19 @@ function set_initial_discrete_controlled_parameters!( (; p) = integrator (; discrete_control) = p - n_conditions = length(discrete_control.condition_value) + n_conditions = sum(length(vec) for vec in discrete_control.condition_value; init = 0) condition_diffs = zeros(Float64, n_conditions) discrete_control_condition(condition_diffs, storage0, integrator.t, integrator) - discrete_control.condition_value .= (condition_diffs .> 0.0) + + # Set the discrete control value (bool) per compound variable + idx_start = 1 + for (compound_variable_idx, vec) in enumerate(discrete_control.condition_value) + l = length(vec) + idx_end = idx_start + l - 1 + discrete_control.condition_value[compound_variable_idx] .= + (condition_diffs[idx_start:idx_end] .> 0) + idx_start += l + end # For every discrete_control node find a condition_idx it listens to for discrete_control_node_id in unique(discrete_control.node_id) @@ -78,7 +87,7 @@ function create_callbacks( saved = SavedResults(saved_flow, saved_vertical_flux, saved_subgrid_level) - n_conditions = length(discrete_control.node_id) + n_conditions = sum(length(vec) for vec in discrete_control.greater_than; init = 0) if n_conditions > 0 discrete_control_cb = VectorContinuousCallback( discrete_control_condition, @@ -183,18 +192,27 @@ Listens for changes in condition truths. function discrete_control_condition(out, u, t, integrator) (; p) = integrator (; discrete_control) = p - - for (i, (listen_node_id, variable, greater_than, look_ahead)) in enumerate( - zip( - discrete_control.listen_node_id, - discrete_control.variable, - discrete_control.greater_than, - discrete_control.look_ahead, - ), + condition_idx = 0 + + # Loop over compound variables + for (listen_node_ids, variables, weights, greater_thans, look_aheads) in zip( + discrete_control.listen_node_id, + discrete_control.variable, + discrete_control.weight, + discrete_control.greater_than, + discrete_control.look_ahead, ) - value = get_value(p, listen_node_id, variable, look_ahead, u, t) - diff = value - greater_than - out[i] = diff + value = 0.0 + for (listen_node_id, variable, weight, look_ahead) in + zip(listen_node_ids, variables, weights, look_aheads) + value += weight * get_value(p, listen_node_id, variable, look_ahead, u, t) + end + # Loop over greater_than values for this compound_variable + for greater_than in greater_thans + condition_idx += 1 + diff = value - greater_than + out[condition_idx] = diff + end end end @@ -252,7 +270,9 @@ function discrete_control_affect_upcrossing!(integrator, condition_idx) (; discrete_control, basin) = p (; variable, condition_value, listen_node_id) = discrete_control - condition_value[condition_idx] = true + compound_variable_idx, greater_than_idx = + get_discrete_control_indices(discrete_control, condition_idx) + condition_value[compound_variable_idx][greater_than_idx] = true control_state_change = discrete_control_affect!(integrator, condition_idx, true) @@ -262,19 +282,24 @@ function discrete_control_affect_upcrossing!(integrator, condition_idx) # only possibly the du. Parameter changes can change the flow on an edge discontinuously, # giving the possibility of logical paradoxes where certain parameter changes immediately # undo the truth state that caused that parameter change. - is_basin = id_index(basin.node_id, discrete_control.listen_node_id[condition_idx])[1] + listen_node_ids = discrete_control.listen_node_id[compound_variable_idx] + is_basin = + length(listen_node_ids) == 1 ? id_index(basin.node_id, only(listen_node_ids))[1] : + false + # NOTE: The above no longer works when listen feature ids can be something other than node ids # I think the more durable option is to give all possible condition types a different variable string, # e.g. basin.level and level_boundary.level - if variable[condition_idx] == "level" && control_state_change && is_basin + if variable[compound_variable_idx][1] == "level" && control_state_change && is_basin # Calling water_balance is expensive, but it is a sure way of getting # du for the basin of this level condition du = zero(u) water_balance!(du, u, p, t) - _, condition_basin_idx = id_index(basin.node_id, listen_node_id[condition_idx]) + _, condition_basin_idx = + id_index(basin.node_id, listen_node_id[compound_variable_idx][1]) if du[condition_basin_idx] < 0.0 - condition_value[condition_idx] = false + condition_value[compound_variable_idx][greater_than_idx] = false discrete_control_affect!(integrator, condition_idx, false) end end @@ -288,7 +313,9 @@ function discrete_control_affect_downcrossing!(integrator, condition_idx) (; discrete_control, basin) = p (; variable, condition_value, listen_node_id) = discrete_control - condition_value[condition_idx] = false + compound_variable_idx, greater_than_idx = + get_discrete_control_indices(discrete_control, condition_idx) + condition_value[compound_variable_idx][greater_than_idx] = false control_state_change = discrete_control_affect!(integrator, condition_idx, false) @@ -298,16 +325,23 @@ function discrete_control_affect_downcrossing!(integrator, condition_idx) # only possibly the du. Parameter changes can change the flow on an edge discontinuously, # giving the possibility of logical paradoxes where certain parameter changes immediately # undo the truth state that caused that parameter change. - if variable[condition_idx] == "level" && control_state_change + compound_variable_idx, greater_than_idx = + get_discrete_control_indices(discrete_control, condition_idx) + listen_node_ids = discrete_control.listen_node_id[compound_variable_idx] + is_basin = + length(listen_node_ids) == 1 ? id_index(basin.node_id, only(listen_node_ids))[1] : + false + + if variable[compound_variable_idx][1] == "level" && control_state_change && is_basin # Calling water_balance is expensive, but it is a sure way of getting # du for the basin of this level condition du = zero(u) water_balance!(du, u, p, t) has_index, condition_basin_idx = - id_index(basin.node_id, listen_node_id[condition_idx]) + id_index(basin.node_id, listen_node_id[compound_variable_idx][1]) if has_index && du[condition_basin_idx] > 0.0 - condition_value[condition_idx] = true + condition_value[compound_variable_idx][greater_than_idx] = true discrete_control_affect!(integrator, condition_idx, true) end end @@ -325,20 +359,34 @@ function discrete_control_affect!( (; discrete_control, graph) = p # Get the discrete_control node that listens to this condition - discrete_control_node_id = discrete_control.node_id[condition_idx] + + compound_variable_idx, _ = get_discrete_control_indices(discrete_control, condition_idx) + discrete_control_node_id = discrete_control.node_id[compound_variable_idx] # Get the indices of all conditions that this control node listens to - condition_ids = discrete_control.node_id .== discrete_control_node_id + where_node_id = searchsorted(discrete_control.node_id, discrete_control_node_id) # Get the truth state for this discrete_control node - truth_values = [ifelse(b, "T", "F") for b in discrete_control.condition_value] - truth_state = join(truth_values[condition_ids], "") + truth_values = cat( + [ + [ifelse(b, "T", "F") for b in discrete_control.condition_value[i]] for + i in where_node_id + ]...; + dims = 1, + ) + truth_state = join(truth_values, "") # Get the truth specific about the latest crossing if !ismissing(upcrossing) - truth_values[condition_idx] = upcrossing ? "U" : "D" + truth_value_idx = + condition_idx - sum( + length(vec) for + vec in discrete_control.condition_value[1:(where_node_id.start - 1)]; + init = 0, + ) + truth_values[truth_value_idx] = upcrossing ? "U" : "D" end - truth_state_crossing_specific = join(truth_values[condition_ids], "") + truth_state_crossing_specific = join(truth_values, "") # What the local control state should be control_state_new = @@ -359,7 +407,7 @@ function discrete_control_affect!( discrete_control.logic_mapping[(discrete_control_node_id, truth_state)] else error( - "Control state specified for neither $truth_state_crossing_specific nor $truth_state for DiscreteControl node $discrete_control_node_id.", + "Control state specified for neither $truth_state_crossing_specific nor $truth_state for $discrete_control_node_id.", ) end diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 16774df21..9ea874d74 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -443,23 +443,28 @@ struct Terminal <: AbstractParameterNode end """ -node_id: node ID of the DiscreteControl node; these are not unique but repeated - by the amount of conditions of this DiscreteControl node -listen_node_id: the ID of the node being condition on -variable: the name of the variable in the condition -greater_than: The threshold value in the condition -condition_value: The current value of each condition +node_id: node ID of the DiscreteControl node per compound variable (can contain repeats) +listen_node_id: the IDs of the nodes being condition on per compound variable +variable: the names of the variables in the condition per compound variable +weight: the weight of the variables in the condition per compound variable +look_ahead: the look ahead of variables in the condition in seconds per compound_variable +greater_than: The threshold values per compound variable +condition_value: The current truth value of each condition per compound_variable per greater_than control_state: Dictionary: node ID => (control state, control state start) logic_mapping: Dictionary: (control node ID, truth state) => control state record: Namedtuple with discrete control information for results """ struct DiscreteControl <: AbstractParameterNode node_id::Vector{NodeID} - listen_node_id::Vector{NodeID} - variable::Vector{String} - look_ahead::Vector{Float64} - greater_than::Vector{Float64} - condition_value::Vector{Bool} + # Definition of compound variables + listen_node_id::Vector{Vector{NodeID}} + variable::Vector{Vector{String}} + weight::Vector{Vector{Float64}} + look_ahead::Vector{Vector{Float64}} + # Definition of conditions (one or more greater_than per compound variable) + greater_than::Vector{Vector{Float64}} + condition_value::Vector{BitVector} + # Definition of logic control_state::Dict{NodeID, Tuple{String, Float64}} logic_mapping::Dict{Tuple{NodeID, String}, String} record::@NamedTuple{ diff --git a/core/src/read.jl b/core/src/read.jl index 20643c943..bd07fd17d 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -542,10 +542,81 @@ function Basin(db::DB, config::Config, chunk_sizes::Vector{Int})::Basin ) end +function parse_variables_and_conditions(compound_variable, condition) + node_id = NodeID[] + listen_node_id = Vector{NodeID}[] + variable = Vector{String}[] + weight = Vector{Float64}[] + look_ahead = Vector{Float64}[] + greater_than = Vector{Float64}[] + condition_value = BitVector[] + errors = false + + # Loop over unique discrete_control node IDs (on which at least one condition is defined) + for id in unique(condition.node_id) + condition_group_id = filter(row -> row.node_id == id, condition) + variable_group_id = filter(row -> row.node_id == id, compound_variable) + # Loop over compound variables for this node ID + for compound_variable_id in unique(condition_group_id.compound_variable_id) + condition_group_variable = filter( + row -> row.compound_variable_id == compound_variable_id, + condition_group_id, + ) + variable_group_variable = filter( + row -> row.compound_variable_id == compound_variable_id, + variable_group_id, + ) + discrete_control_id = NodeID(NodeType.DiscreteControl, id) + if isempty(variable_group_variable) + errors = true + @error "compound_variable_id $compound_variable_id for $discrete_control_id in condition table but not in variable table" + else + push!(node_id, discrete_control_id) + push!( + listen_node_id, + NodeID.( + variable_group_variable.listen_node_type, + variable_group_variable.listen_node_id, + ), + ) + push!(variable, variable_group_variable.variable) + push!(weight, coalesce.(variable_group_variable.weight, 1.0)) + push!(look_ahead, coalesce.(variable_group_variable.look_ahead, 0.0)) + push!(greater_than, condition_group_variable.greater_than) + push!( + condition_value, + BitVector(zeros(length(condition_group_variable.greater_than))), + ) + end + end + end + return node_id, + listen_node_id, + variable, + weight, + look_ahead, + greater_than, + condition_value, + !errors +end + function DiscreteControl(db::DB, config::Config)::DiscreteControl condition = load_structvector(db, config, DiscreteControlConditionV1) + compound_variable = load_structvector(db, config, DiscreteControlVariableV1) + + node_id, + listen_node_id, + variable, + weight, + look_ahead, + greater_than, + condition_value, + valid = parse_variables_and_conditions(compound_variable, condition) + + if !valid + error("Problems encountered when parsing DiscreteControl variables and conditions.") + end - condition_value = fill(false, length(condition.node_id)) control_state::Dict{NodeID, Tuple{String, Float64}} = Dict() rows = execute(db, "SELECT from_node_id, edge_type FROM Edge ORDER BY fid") @@ -557,7 +628,6 @@ function DiscreteControl(db::DB, config::Config)::DiscreteControl end logic = load_structvector(db, config, DiscreteControlLogicV1) - logic_mapping = Dict{Tuple{NodeID, String}, String}() for (node_id, truth_state, control_state_) in @@ -567,7 +637,6 @@ function DiscreteControl(db::DB, config::Config)::DiscreteControl end logic_mapping = expand_logic_mapping(logic_mapping) - look_ahead = coalesce.(condition.look_ahead, 0.0) record = ( time = Float64[], @@ -577,11 +646,12 @@ function DiscreteControl(db::DB, config::Config)::DiscreteControl ) return DiscreteControl( - NodeID.(NodeType.DiscreteControl, condition.node_id), # Not unique - NodeID.(condition.listen_node_type, condition.listen_node_id), - condition.variable, + node_id, # Not unique + listen_node_id, + variable, + weight, look_ahead, - condition.greater_than, + greater_than, condition_value, control_state, logic_mapping, diff --git a/core/src/schema.jl b/core/src/schema.jl index d59190c3b..0875e48c5 100644 --- a/core/src/schema.jl +++ b/core/src/schema.jl @@ -1,5 +1,6 @@ # These schemas define the name of database tables and the configuration file structure # The identifier is parsed as ribasim.nodetype.kind, no capitals or underscores are allowed. +@schema "ribasim.discretecontrol.variable" DiscreteControlVariable @schema "ribasim.discretecontrol.condition" DiscreteControlCondition @schema "ribasim.discretecontrol.logic" DiscreteControlLogic @schema "ribasim.basin.static" BasinStatic @@ -183,15 +184,22 @@ end node_id::Int32 end -@version DiscreteControlConditionV1 begin +@version DiscreteControlVariableV1 begin node_id::Int32 + compound_variable_id::Int32 listen_node_type::String listen_node_id::Int32 variable::String - greater_than::Float64 + weight::Union{Missing, Float64} look_ahead::Union{Missing, Float64} end +@version DiscreteControlConditionV1 begin + node_id::Int32 + compound_variable_id::Int32 + greater_than::Float64 +end + @version DiscreteControlLogicV1 begin node_id::Int32 truth_state::String diff --git a/core/src/util.jl b/core/src/util.jl index 46a3643b9..57482cca1 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -709,3 +709,19 @@ function get_influx(basin::Basin, basin_idx::Int)::Float64 return precipitation[basin_idx] - evaporation[basin_idx] + drainage[basin_idx] - infiltration[basin_idx] end + +function get_discrete_control_indices(discrete_control::DiscreteControl, condition_idx::Int) + (; greater_than) = discrete_control + condition_idx_now = 1 + + for (compound_variable_idx, vec) in enumerate(greater_than) + l = length(vec) + + if condition_idx_now + l > condition_idx + greater_than_idx = condition_idx - condition_idx_now + 1 + return compound_variable_idx, greater_than_idx + end + + condition_idx_now += l + end +end diff --git a/core/src/validation.jl b/core/src/validation.jl index 72fcb55ff..ae039a561 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -102,6 +102,9 @@ sort_by_id_state_level(row) = (row.node_id, row.control_state, row.level) sort_by_priority(row) = (row.node_id, row.priority) sort_by_priority_time(row) = (row.node_id, row.priority, row.time) sort_by_subgrid_level(row) = (row.subgrid_id, row.basin_level) +sort_by_variable(row) = + (row.node_id, row.listen_node_type, row.listen_node_id, row.variable) +sort_by_condition(row) = (row.node_id, row.compound_variable_id, row.greater_than) # get the right sort by function given the Schema, with sort_by_id as the default sort_by_function(table::StructVector{<:Legolas.AbstractRecord}) = sort_by_id @@ -110,6 +113,8 @@ sort_by_function(table::StructVector{BasinProfileV1}) = sort_by_id_level sort_by_function(table::StructVector{UserDemandStaticV1}) = sort_by_priority sort_by_function(table::StructVector{UserDemandTimeV1}) = sort_by_priority_time sort_by_function(table::StructVector{BasinSubgridV1}) = sort_by_subgrid_level +sort_by_function(table::StructVector{DiscreteControlVariableV1}) = sort_by_variable +sort_by_function(table::StructVector{DiscreteControlConditionV1}) = sort_by_condition const TimeSchemas = Union{ BasinTimeV1, @@ -502,7 +507,8 @@ Check: """ function valid_discrete_control(p::Parameters, config::Config)::Bool (; discrete_control, graph) = p - (; node_id, logic_mapping, look_ahead, variable, listen_node_id) = discrete_control + (; node_id, logic_mapping, look_ahead, variable, listen_node_id, greater_than) = + discrete_control t_end = seconds_since(config.endtime, config.starttime) errors = false @@ -515,7 +521,7 @@ function valid_discrete_control(p::Parameters, config::Config)::Bool truth_states_wrong_length = String[] # The number of conditions of this DiscreteControl node - n_conditions = length(searchsorted(node_id, id)) + n_conditions = sum(length(greater_than[i]) for i in searchsorted(node_id, id)) for (key, control_state) in logic_mapping id_, truth_state = key @@ -563,29 +569,30 @@ function valid_discrete_control(p::Parameters, config::Config)::Bool end end end - for (Δt, var, node_id) in zip(look_ahead, variable, listen_node_id) - if !iszero(Δt) - node_type = node_id.type - # TODO: If more transient listen variables must be supported, this validation must be more specific - # (e.g. for some node some variables are transient, some not). - if node_type ∉ [NodeType.FlowBoundary, NodeType.LevelBoundary] - errors = true - @error "Look ahead supplied for non-timeseries listen variable '$var' from listen node $node_id." - else - if Δt < 0 + for (look_aheads, variables, listen_node_ids) in + zip(look_ahead, variable, listen_node_id) + for (Δt, var, node_id) in zip(look_aheads, variables, listen_node_ids) + if !iszero(Δt) + node_type = node_id.type + if node_type ∉ [NodeType.FlowBoundary, NodeType.LevelBoundary] errors = true - @error "Negative look ahead supplied for listen variable '$var' from listen node $node_id." + @error "Look ahead supplied for non-timeseries listen variable '$var' from listen node $node_id." else - node = getfield(p, graph[node_id].type) - idx = if node_type == NodeType.Basin - id_index(node.node_id, node_id) - else - searchsortedfirst(node.node_id, node_id) - end - interpolation = getfield(node, Symbol(var))[idx] - if t_end + Δt > interpolation.t[end] + if Δt < 0 errors = true - @error "Look ahead for listen variable '$var' from listen node $node_id goes past timeseries end during simulation." + @error "Negative look ahead supplied for listen variable '$var' from listen node $node_id." + else + node = getfield(p, graph[node_id].type) + idx = if node_type == NodeType.Basin + id_index(node.node_id, node_id) + else + searchsortedfirst(node.node_id, node_id) + end + interpolation = getfield(node, Symbol(var))[idx] + if t_end + Δt > interpolation.t[end] + errors = true + @error "Look ahead for listen variable '$var' from listen node $node_id goes past timeseries end during simulation." + end end end end diff --git a/core/test/control_test.jl b/core/test/control_test.jl index fdb2e80f8..21bcc653c 100644 --- a/core/test/control_test.jl +++ b/core/test/control_test.jl @@ -37,11 +37,11 @@ # Control times t_1 = discrete_control.record.time[3] t_1_index = findfirst(>=(t_1), t) - @test level[1, t_1_index] <= discrete_control.greater_than[1] + @test level[1, t_1_index] <= discrete_control.greater_than[1][1] t_2 = discrete_control.record.time[4] t_2_index = findfirst(>=(t_2), t) - @test level[2, t_2_index] >= discrete_control.greater_than[2] + @test level[2, t_2_index] >= discrete_control.greater_than[2][1] flow = get_tmp(graph[].flow, 0) @test all(iszero, flow) @@ -54,13 +54,13 @@ end p = model.integrator.p (; discrete_control, flow_boundary) = p - Δt = discrete_control.look_ahead[1] + Δt = discrete_control.look_ahead[1][1] t = Ribasim.tsaves(model) t_control = discrete_control.record.time[2] t_control_index = searchsortedfirst(t, t_control) - greater_than = discrete_control.greater_than[1] + greater_than = discrete_control.greater_than[1][1] flow_t_control = flow_boundary.flow_rate[1](t_control) flow_t_control_ahead = flow_boundary.flow_rate[1](t_control + Δt) @@ -78,13 +78,13 @@ end p = model.integrator.p (; discrete_control, level_boundary) = p - Δt = discrete_control.look_ahead[1] + Δt = discrete_control.look_ahead[1][1] t = Ribasim.tsaves(model) t_control = discrete_control.record.time[2] t_control_index = searchsortedfirst(t, t_control) - greater_than = discrete_control.greater_than[1] + greater_than = discrete_control.greater_than[1][1] level_t_control = level_boundary.level[1](t_control) level_t_control_ahead = level_boundary.level[1](t_control + Δt) @@ -167,8 +167,8 @@ end t_in = discrete_control.record.time[3] t_none_2 = discrete_control.record.time[4] - level_min = greater_than[1] - setpoint = greater_than[2] + level_min = greater_than[1][1] + setpoint = greater_than[1][2] t_1_none_index = findfirst(>=(t_none_1), t) t_in_index = findfirst(>=(t_in), t) @@ -206,3 +206,23 @@ end @test isapprox(level[t_idx_target_jump], target_high, atol = 1e-1) @test isapprox(level[end], target_low, atol = 1e-1) end + +@testitem "Compound condition" begin + using Ribasim: NodeID + + toml_path = normpath( + @__DIR__, + "../../generated_testmodels/compound_variable_condition/ribasim.toml", + ) + @test ispath(toml_path) + model = Ribasim.run(toml_path) + (; discrete_control) = model.integrator.p + (; listen_node_id, variable, weight, record) = discrete_control + + @test listen_node_id == [[NodeID(:FlowBoundary, 2), NodeID(:FlowBoundary, 3)]] + @test variable == [["flow_rate", "flow_rate"]] + @test weight == [[0.5, 0.5]] + @test record.time ≈ [0.0, model.integrator.sol.t[end] / 2] + @test record.truth_state == ["F", "T"] + @test record.control_state == ["Off", "On"] +end diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index 3613fa409..88374f8ef 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -204,7 +204,8 @@ name it must have in the database if it is stored there. - 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) - - `DisceteControl / condition`: Conditions of the form 'the level in the basin with node id `n` is bigger than 2.0 m' + - `DiscreteControl / variable` : Linear combinations of variables, e.g. averages or differences of levels + - `DisceteControl / condition`: Conditions of the form 'the compound variable with ID `n` is bigger than 2.0' - `DisceteControl / logic`: Translates the truth value of a set of conditions to parameter values for a controlled node - PidControl: Controls the level in a basin by continuously controlling the flow rate of a connected pump or outlet. See also [Wikipedia](https://en.wikipedia.org/wiki/PID_controller) and [PID controller in equations](equations.qmd#sec-PID). - `PidControl / static`: The proportional, integral, and derivative parameters, the target value and which basin should be controlled @@ -625,18 +626,35 @@ node_id | Int32 | - | sorted DiscreteControl is implemented based on [VectorContinuousCallback](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#VectorContinuousCallback). +## DiscreteControl / 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. + +column | type | unit | restriction +-------------------- | -------- | ------- | ----------- +node_id | Int32 | - | sorted +compound_variable_id | Int32 | - | sorted per node_id +listen_node_type | String | - | known node type +listen_node_id | Int32 | - | sorted per node_id +variable | String | - | must be "level" or "flow_rate", sorted per listen_node_id +weight | Float64 | - | (optional, default 1.0) +look_ahead | Float64 | $s$ | Only on transient boundary conditions, non-negative (optional, default 0.0). + ## DiscreteControl / condition -The condition schema defines conditions of the form 'the discrete_control node with this node id listens to whether the given variable of the node with the given listen feature id is grater than the given value'. If the condition variable comes from a time-series, a look ahead $\Delta t$ can be supplied. +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 `greater_than`'. Multiple conditions with different `greater_than` values can be defined on the same compound_variable. -column | type | unit | restriction ------------------ | -------- | ------- | ----------- -node_id | Int32 | - | sorted -listen_node_type | String | - | known node type -listen_node_id | Int32 | - | sorted per node_id -variable | String | - | must be "level" or "flow_rate", sorted per listen_node_id -greater_than | Float64 | various | sorted per variable -look_ahead | Float64 | $s$ | Only on transient boundary conditions, non-negative (optional, default 0) +column | type | unit | restriction +-------------------- | -------- | ------- | ----------- +node_id | Int32 | - | sorted +compound_variable_id | Int32 | - | - +greater_than | Float64 | various | sorted per variable ## DiscreteControl / logic diff --git a/docs/python/examples.ipynb b/docs/python/examples.ipynb index 5a76e3ae7..32bbb556b 100644 --- a/docs/python/examples.ipynb +++ b/docs/python/examples.ipynb @@ -492,11 +492,14 @@ "model.discrete_control.add(\n", " Node(7, Point(1.0, 0.0)),\n", " [\n", + " discrete_control.Variable(\n", + " listen_node_id=[1],\n", + " listen_node_type=\"Basin\",\n", + " variable=\"level\",\n", + " compound_variable_id=1,\n", + " ),\n", " discrete_control.Condition(\n", - " listen_node_id=[1, 1, 1],\n", - " listen_node_type=[\"Basin\", \"Basin\", \"Basin\"],\n", - " variable=[\"level\", \"level\", \"level\"],\n", - " greater_than=[5.0, 10.0, 15.0],\n", + " greater_than=[5.0, 10.0, 15.0], compound_variable_id=1\n", " ),\n", " discrete_control.Logic(\n", " truth_state=[\"FFF\", \"U**\", \"T*F\", \"**D\", \"TTT\"],\n", @@ -691,7 +694,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now run the model with `level_setpoint_with_minmax/ribasim.toml`.\n", + "Now run the model with `ribasim level_setpoint_with_minmax/ribasim.toml`.\n", "After running the model, read back the results:\n" ] }, @@ -1183,10 +1186,14 @@ "model.discrete_control.add(\n", " Node(11, Point(4.5, 0.25), subnetwork_id=1),\n", " [\n", - " discrete_control.Condition(\n", + " discrete_control.Variable(\n", + " compound_variable_id=1,\n", " listen_node_id=[5],\n", " listen_node_type=[\"Basin\"],\n", " variable=[\"level\"],\n", + " ),\n", + " discrete_control.Condition(\n", + " compound_variable_id=1,\n", " greater_than=[0.52],\n", " ),\n", " discrete_control.Logic(\n", diff --git a/python/ribasim/ribasim/config.py b/python/ribasim/ribasim/config.py index f9605111d..30408fe05 100644 --- a/python/ribasim/ribasim/config.py +++ b/python/ribasim/ribasim/config.py @@ -23,6 +23,7 @@ BasinTimeSchema, DiscreteControlConditionSchema, DiscreteControlLogicSchema, + DiscreteControlVariableSchema, FlowBoundaryStaticSchema, FlowBoundaryTimeSchema, FlowDemandStaticSchema, @@ -287,10 +288,25 @@ class ManningResistance(MultiNodeModel): class DiscreteControl(MultiNodeModel): + variable: TableModel[DiscreteControlVariableSchema] = Field( + default_factory=TableModel[DiscreteControlVariableSchema], + json_schema_extra={ + "sort_keys": [ + "node_id", + "listen_node_type", + "listen_node_id", + "variable", + ] + }, + ) condition: TableModel[DiscreteControlConditionSchema] = Field( default_factory=TableModel[DiscreteControlConditionSchema], json_schema_extra={ - "sort_keys": ["node_id", "listen_node_id", "variable", "greater_than"] + "sort_keys": [ + "node_id", + "compound_variable_id", + "greater_than", + ] }, ) logic: TableModel[DiscreteControlLogicSchema] = Field( diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py index 895cf7220..2b4f195ef 100644 --- a/python/ribasim/ribasim/model.py +++ b/python/ribasim/ribasim/model.py @@ -321,9 +321,9 @@ def plot_control_listen(self, ax): df_listen_edge = pd.concat([df_listen_edge, to_add]) # Listen edges from DiscreteControl - condition = self.discrete_control.condition.df - if condition is not None: - to_add = condition[ + df_variable = self.discrete_control.variable.df + if df_variable is not None: + to_add = df_variable[ ["node_id", "listen_node_id", "listen_node_type"] ].drop_duplicates() to_add.columns = ["control_node_id", "listen_node_id", "listen_node_type"] diff --git a/python/ribasim/ribasim/nodes/discrete_control.py b/python/ribasim/ribasim/nodes/discrete_control.py index eeb0b28ca..df7b9f42b 100644 --- a/python/ribasim/ribasim/nodes/discrete_control.py +++ b/python/ribasim/ribasim/nodes/discrete_control.py @@ -1,9 +1,18 @@ from pandas import DataFrame from ribasim.input_base import TableModel -from ribasim.schemas import DiscreteControlConditionSchema, DiscreteControlLogicSchema +from ribasim.schemas import ( + DiscreteControlConditionSchema, + DiscreteControlLogicSchema, + DiscreteControlVariableSchema, +) -__all__ = ["Condition", "Logic"] +__all__ = ["Condition", "Logic", "Variable"] + + +class Variable(TableModel[DiscreteControlVariableSchema]): + def __init__(self, **kwargs): + super().__init__(df=DataFrame(dict(**kwargs))) class Condition(TableModel[DiscreteControlConditionSchema]): diff --git a/python/ribasim/ribasim/schemas.py b/python/ribasim/ribasim/schemas.py index 0cf033671..2a8b5eab5 100644 --- a/python/ribasim/ribasim/schemas.py +++ b/python/ribasim/ribasim/schemas.py @@ -50,11 +50,8 @@ class BasinTimeSchema(_BaseSchema): class DiscreteControlConditionSchema(_BaseSchema): node_id: Series[Int32] = pa.Field(nullable=False, default=0) - listen_node_type: Series[str] = pa.Field(nullable=False) - listen_node_id: Series[Int32] = pa.Field(nullable=False, default=0) - variable: Series[str] = pa.Field(nullable=False) + compound_variable_id: Series[Int32] = pa.Field(nullable=False, default=0) greater_than: Series[float] = pa.Field(nullable=False) - look_ahead: Series[float] = pa.Field(nullable=True) class DiscreteControlLogicSchema(_BaseSchema): @@ -63,6 +60,16 @@ class DiscreteControlLogicSchema(_BaseSchema): control_state: Series[str] = pa.Field(nullable=False) +class DiscreteControlVariableSchema(_BaseSchema): + node_id: Series[Int32] = pa.Field(nullable=False, default=0) + compound_variable_id: Series[Int32] = pa.Field(nullable=False, default=0) + listen_node_type: Series[str] = pa.Field(nullable=False) + listen_node_id: Series[Int32] = pa.Field(nullable=False, default=0) + variable: Series[str] = pa.Field(nullable=False) + weight: Series[float] = pa.Field(nullable=True) + look_ahead: Series[float] = pa.Field(nullable=True) + + class FlowBoundaryStaticSchema(_BaseSchema): node_id: Series[Int32] = pa.Field(nullable=False, default=0) active: Series[pa.BOOL] = pa.Field(nullable=True) diff --git a/python/ribasim/tests/test_io.py b/python/ribasim/tests/test_io.py index bd6fcc7a8..e6507d77e 100644 --- a/python/ribasim/tests/test_io.py +++ b/python/ribasim/tests/test_io.py @@ -102,8 +102,7 @@ def test_sort(level_setpoint_with_minmax, tmp_path): assert table.df.iloc[0]["greater_than"] == 15.0 assert table._sort_keys == [ "node_id", - "listen_node_id", - "variable", + "compound_variable_id", "greater_than", ] table.sort() diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index ceeb1fad2..4f9ce86c3 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -28,6 +28,7 @@ ) from ribasim_testmodels.bucket import bucket_model, leaky_bucket_model from ribasim_testmodels.discrete_control import ( + compound_variable_condition_model, flow_condition_model, level_boundary_condition_model, level_setpoint_with_minmax_model, @@ -64,6 +65,7 @@ "basic_model", "basic_transient_model", "bucket_model", + "compound_variable_condition_model", "discrete_control_of_pid_control_model", "dutch_waterways_model", "flow_boundary_time_model", diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 1e102e8c0..51980abc5 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -406,11 +406,15 @@ def fractional_flow_subnetwork_model() -> Model: model.discrete_control.add( Node(10, Point(-1, 2), subnetwork_id=2), [ - discrete_control.Condition( + discrete_control.Variable( listen_node_type="FlowBoundary", listen_node_id=[1], variable="flow_rate", - greater_than=3e-3, + compound_variable_id=1, + ), + discrete_control.Condition( + greater_than=[3e-3], + compound_variable_id=1, ), discrete_control.Logic(truth_state=["F", "T"], control_state=["A", "B"]), ], @@ -500,11 +504,15 @@ def allocation_example_model() -> Model: model.discrete_control.add( Node(11, Point(4.5, 0.25), subnetwork_id=2), [ - discrete_control.Condition( + discrete_control.Variable( listen_node_type="Basin", listen_node_id=[5], variable="level", - greater_than=0.52, + compound_variable_id=1, + ), + discrete_control.Condition( + greater_than=[0.52], + compound_variable_id=1, ), discrete_control.Logic( truth_state=["T", "F"], control_state=["divert", "close"] @@ -679,11 +687,15 @@ def main_network_with_subnetworks_model() -> Model: model.discrete_control.add( Node(33, Point(13, 5), subnetwork_id=5), [ - discrete_control.Condition( + discrete_control.Variable( listen_node_type="Basin", listen_node_id=[25], variable="level", - greater_than=0.003, + compound_variable_id=1, + ), + discrete_control.Condition( + greater_than=[0.003], + compound_variable_id=1, ), discrete_control.Logic(truth_state=["F", "T"], control_state=["A", "B"]), ], diff --git a/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py b/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py index aae7674c9..650e73876 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py +++ b/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py @@ -54,11 +54,15 @@ def pump_discrete_control_model() -> Model: model.discrete_control.add( Node(5, Point(1, 1)), [ - discrete_control.Condition( + discrete_control.Variable( listen_node_type="Basin", listen_node_id=[1, 3], variable="level", + compound_variable_id=[1, 2], + ), + discrete_control.Condition( greater_than=[0.8, 0.4], + compound_variable_id=[1, 2], ), discrete_control.Logic( truth_state=["FF", "TF", "FT", "TT"], @@ -69,11 +73,15 @@ def pump_discrete_control_model() -> Model: model.discrete_control.add( Node(6, Point(2, -1)), [ - discrete_control.Condition( + discrete_control.Variable( listen_node_type="Basin", - listen_node_id=3, + listen_node_id=[3], variable="level", + compound_variable_id=1, + ), + discrete_control.Condition( greater_than=[0.45], + compound_variable_id=1, ), discrete_control.Logic( truth_state=["T", "F"], @@ -140,12 +148,16 @@ def flow_condition_model() -> Model: model.discrete_control.add( Node(5, Point(1, 1)), [ - discrete_control.Condition( + discrete_control.Variable( listen_node_type="FlowBoundary", - listen_node_id=1, + listen_node_id=[1], variable="flow_rate", - greater_than=[20 / (86400)], look_ahead=60 * 86400, + compound_variable_id=1, + ), + discrete_control.Condition( + greater_than=[20 / (86400)], + compound_variable_id=1, ), discrete_control.Logic(truth_state=["T", "F"], control_state=["off", "on"]), ], @@ -203,12 +215,16 @@ def level_boundary_condition_model() -> Model: model.discrete_control.add( Node(6, Point(1.5, 1)), [ - discrete_control.Condition( + discrete_control.Variable( listen_node_type="LevelBoundary", listen_node_id=[1], variable="level", - greater_than=6.0, look_ahead=60 * 86400, + compound_variable_id=1, + ), + discrete_control.Condition( + greater_than=[6.0], + compound_variable_id=1, ), discrete_control.Logic(truth_state=["T", "F"], control_state=["on", "off"]), ], @@ -275,11 +291,15 @@ def tabulated_rating_curve_control_model() -> Model: model.discrete_control.add( Node(4, Point(1, 1)), [ - discrete_control.Condition( + discrete_control.Variable( listen_node_type="Basin", listen_node_id=[1], variable="level", - greater_than=0.5, + compound_variable_id=1, + ), + discrete_control.Condition( + greater_than=[0.5], + compound_variable_id=1, ), discrete_control.Logic( truth_state=["T", "F"], control_state=["low", "high"] @@ -342,12 +362,16 @@ def level_setpoint_with_minmax_model() -> Model: model.discrete_control.add( Node(7, Point(1, 0)), [ - discrete_control.Condition( + discrete_control.Variable( listen_node_type="Basin", - listen_node_id=1, + listen_node_id=[1], variable="level", + compound_variable_id=1, + ), + discrete_control.Condition( # min, setpoint, max greater_than=[5.0, 10.0, 15.0], + compound_variable_id=1, ), discrete_control.Logic( truth_state=["FFF", "U**", "T*F", "**D", "TTT"], @@ -390,3 +414,61 @@ def level_setpoint_with_minmax_model() -> Model: ) return model + + +def compound_variable_condition_model() -> Model: + """ + Set up a minimal model containing a condition on a compound variable + for discrete control. + """ + + model = Model( + starttime="2020-01-01", + endtime="2021-01-01", + crs="EPSG:28992", + ) + + model.basin.add( + Node(1, Point(1, 0)), + [ + basin.Profile(area=1000.0, level=[0.0, 1.0]), + basin.State(level=[1.0]), + ], + ) + model.flow_boundary.add( + Node(2, Point(0, 0)), [flow_boundary.Static(flow_rate=[0.0])] + ) + model.flow_boundary.add( + Node(3, Point(0, 1)), + [flow_boundary.Time(time=["2020-01-01", "2021-01-01"], flow_rate=[0.0, 2.0])], + ) + model.pump.add( + Node(4, Point(2, 0)), + [pump.Static(control_state=["Off", "On"], flow_rate=[0.0, 1.0])], + ) + model.terminal.add(Node(5, Point(3, 0))) + model.discrete_control.add( + Node(6, Point(1, 1)), + [ + discrete_control.Variable( + listen_node_type="FlowBoundary", + listen_node_id=[2, 3], + variable="flow_rate", + weight=0.5, + compound_variable_id=1, + ), + discrete_control.Condition( + greater_than=[0.5], + compound_variable_id=1, + ), + discrete_control.Logic(truth_state=["T", "F"], control_state=["On", "Off"]), + ], + ) + + model.edge.add(model.flow_boundary[2], model.basin[1]) + model.edge.add(model.flow_boundary[3], model.basin[1]) + model.edge.add(model.basin[1], model.pump[4]) + model.edge.add(model.pump[4], model.terminal[5]) + model.edge.add(model.discrete_control[6], model.pump[4]) + + return model diff --git a/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py b/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py index fd990957b..9fb315a7f 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py +++ b/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py @@ -121,11 +121,15 @@ def dutch_waterways_model() -> Model: model.discrete_control.add( Node(17, Point(183612, 441833), name="Controller Driel"), [ - discrete_control.Condition( + discrete_control.Variable( listen_node_type="FlowBoundary", - listen_node_id=1, + listen_node_id=[1], variable="flow_rate", + compound_variable_id=1, + ), + discrete_control.Condition( greater_than=[250, 275, 750, 800], + compound_variable_id=1, ), discrete_control.Logic( truth_state=["FFFF", "U***", "T**F", "***D", "TTTT"], diff --git a/python/ribasim_testmodels/ribasim_testmodels/invalid.py b/python/ribasim_testmodels/ribasim_testmodels/invalid.py index 906a7b485..463603bf6 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/invalid.py +++ b/python/ribasim_testmodels/ribasim_testmodels/invalid.py @@ -148,15 +148,19 @@ def invalid_discrete_control_model() -> Model: model.discrete_control.add( Node(5, Point(1, 1)), [ - discrete_control.Condition( + discrete_control.Variable( listen_node_type=["Basin", "FlowBoundary", "FlowBoundary"], listen_node_id=[1, 4, 4], variable=["level", "flow_rate", "flow_rate"], - greater_than=[0.5, 1.5, 1.5], # Invalid: look_ahead can only be specified for timeseries variables. # Invalid: this look_ahead will go past the provided timeseries during simulation. # Invalid: look_ahead must be non-negative. look_ahead=[100.0, 40 * 24 * 60 * 60, -10.0], + compound_variable_id=[1, 2, 3], + ), + discrete_control.Condition( + greater_than=[0.5, 1.5, 1.5], + compound_variable_id=[1, 2, 3], ), # Invalid: DiscreteControl node #4 has 2 conditions so # truth states have to be of length 2 diff --git a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py index d8c591c09..4a543acc3 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py +++ b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py @@ -129,11 +129,15 @@ def discrete_control_of_pid_control_model() -> Model: model.discrete_control.add( Node(7, Point(0, 1)), [ - discrete_control.Condition( + discrete_control.Variable( listen_node_type="LevelBoundary", listen_node_id=[1], variable="level", - greater_than=5.0, + compound_variable_id=1, + ), + discrete_control.Condition( + greater_than=[5.0], + compound_variable_id=1, ), discrete_control.Logic( truth_state=["T", "F"], control_state=["target_high", "target_low"] diff --git a/ribasim_qgis/core/nodes.py b/ribasim_qgis/core/nodes.py index 95edb13bd..97682e780 100644 --- a/ribasim_qgis/core/nodes.py +++ b/ribasim_qgis/core/nodes.py @@ -691,10 +691,10 @@ def attributes(cls) -> list[QgsField]: ] -class DiscreteControlCondition(Input): +class DiscreteControlVariable(Input): @classmethod def input_type(cls) -> str: - return "DiscreteControl / condition" + return "DiscreteControl / variable" @classmethod def geometry_type(cls) -> str: @@ -704,9 +704,29 @@ def geometry_type(cls) -> str: def attributes(cls) -> list[QgsField]: return [ QgsField("node_id", QVariant.Int), + QgsField("compound_variable_id", QVariant.Int), QgsField("listen_node_type", QVariant.String), QgsField("listen_node_id", QVariant.Int), QgsField("variable", QVariant.String), + QgsField("weight", QVariant.Double), + QgsField("look_ahead", QVariant.Double), + ] + + +class DiscreteControlCondition(Input): + @classmethod + def input_type(cls) -> str: + return "DiscreteControl / condition" + + @classmethod + def geometry_type(cls) -> str: + return "No Geometry" + + @classmethod + def attributes(cls) -> list[QgsField]: + return [ + QgsField("node_id", QVariant.Int), + QgsField("compound_variable_id", QVariant.Int), QgsField("greater_than", QVariant.Double), ]