diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index a6ec614c9..c862406ef 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -1009,9 +1009,38 @@ function optimize_demand_priority!( demand_priorities_all[demand_priority_idx], optimization_type, ) + + if optimization_type == OptimizationType.allocate + apply_control_from_allocation!(p, allocation_model) + end return nothing end +function apply_control_from_allocation!(p::Parameters, allocation_model::AllocationModel) + (; pump, outlet, graph) = p + + apply_control_from_allocation!(pump, allocation_model, graph) + apply_control_from_allocation!(outlet, allocation_model, graph) +end + +function apply_control_from_allocation!( + node::Union{Pump, Outlet}, + allocation_model::AllocationModel, + graph::MetaGraph, +) + (; node_id, control_type, inflow_link) = node + (; flow, subnetwork_id) = allocation_model + cache = node.flow_rate_cache[Float64[]] # TODO: Make non-allocating + + for (id, c_type, link_in) in zip(node_id, control_type, inflow_link) + in_subnetwork = (graph[id].subnetwork_id == subnetwork_id) + allocation_controlled = (c_type == ControlType.Allocation) + if in_subnetwork && allocation_controlled + cache[id.idx] = flow[link_in.link] + end + end +end + """ Set the initial capacities and demands which are reduced by usage. """ diff --git a/core/src/callback.jl b/core/src/callback.jl index 75cd5e87c..dff1753bd 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -651,7 +651,7 @@ end function apply_parameter_update!(parameter_update)::Nothing (; name, value, ref) = parameter_update - # Ignore this parameter update of the associated node does + # Ignore this parameter update if the associated node does # not have an 'active' field if name == :active && ref.i == 0 return nothing diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 63b196639..6b3c8a6f1 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -16,7 +16,7 @@ const SolverStats = @NamedTuple{ # LinkType.flow and NodeType.FlowBoundary @enumx LinkType flow control none @eval @enumx NodeType $(config.nodetypes...) -@enumx ContinuousControlType None Continuous PID +@enumx ControlType None Continuous PID Allocation @enumx Substance Continuity = 1 Initial = 2 LevelBoundary = 3 FlowBoundary = 4 UserDemand = 5 Drainage = 6 Precipitation = 7 Base.to_index(id::Substance.T) = Int(id) # used to index into concentration matrices @@ -650,7 +650,7 @@ max_flow_rate: The maximum flow rate of the pump min_upstream_level: The upstream level below which the Pump flow goes to zero max_downstream_level: The downstream level above which the Pump flow goes to zero control_mapping: dictionary from (node_id, control_state) to target flow rate -continuous_control_type: one of None, ContinuousControl, PidControl +control_type: one of None, ContinuousControl, PidControl, Allocation """ @kwdef struct Pump <: AbstractParameterNode node_id::Vector{NodeID} @@ -664,8 +664,7 @@ continuous_control_type: one of None, ContinuousControl, PidControl min_upstream_level::Vector{ScalarInterpolation} = ScalarInterpolation[] max_downstream_level::Vector{ScalarInterpolation} = ScalarInterpolation[] control_mapping::Dict{Tuple{NodeID, String}, ControlStateUpdate} - continuous_control_type::Vector{ContinuousControlType.T} = - fill(ContinuousControlType.None, length(node_id)) + control_type::Vector{ControlType.T} = fill(ControlType.None, length(node_id)) function Pump( node_id, @@ -679,7 +678,7 @@ continuous_control_type: one of None, ContinuousControl, PidControl min_upstream_level, max_downstream_level, control_mapping, - continuous_control_type, + control_type, ) if valid_flow_rates(node_id, flow_rate[Float64[]], control_mapping) return new( @@ -694,7 +693,7 @@ continuous_control_type: one of None, ContinuousControl, PidControl min_upstream_level, max_downstream_level, control_mapping, - continuous_control_type, + control_type, ) else error("Invalid Pump flow rate(s).") @@ -716,7 +715,7 @@ max_flow_rate: The maximum flow rate of the outlet min_upstream_level: The upstream level below which the Outlet flow goes to zero max_downstream_level: The downstream level above which the Outlet flow goes to zero control_mapping: dictionary from (node_id, control_state) to target flow rate -continuous_control_type: one of None, ContinuousControl, PidControl +control_type: one of None, ContinuousControl, PidControl, Allocation """ @kwdef struct Outlet <: AbstractParameterNode node_id::Vector{NodeID} @@ -730,8 +729,7 @@ continuous_control_type: one of None, ContinuousControl, PidControl min_upstream_level::Vector{ScalarInterpolation} = ScalarInterpolation[] max_downstream_level::Vector{ScalarInterpolation} = ScalarInterpolation[] control_mapping::Dict{Tuple{NodeID, String}, ControlStateUpdate} = Dict() - continuous_control_type::Vector{ContinuousControlType.T} = - fill(ContinuousControlType.None, length(node_id)) + control_type::Vector{ControlType.T} = fill(ControlType.None, length(node_id)) function Outlet( node_id, @@ -745,7 +743,7 @@ continuous_control_type: one of None, ContinuousControl, PidControl min_upstream_level, max_downstream_level, control_mapping, - continuous_control_type, + control_type, ) if valid_flow_rates(node_id, flow_rate[Float64[]], control_mapping) return new( @@ -760,7 +758,7 @@ continuous_control_type: one of None, ContinuousControl, PidControl min_upstream_level, max_downstream_level, control_mapping, - continuous_control_type, + control_type, ) else error("Invalid Outlet flow rate(s).") diff --git a/core/src/read.jl b/core/src/read.jl index f5ccaa575..b584e595e 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -742,7 +742,10 @@ function Pump(db::DB, config::Config, graph::MetaGraph)::Pump flow_rate_cache = cache(length(node_id)) flow_rate_cache[Float64[]] .= [itp(0) for itp in parsed_parameters.flow_rate] - return Pump(; + control_type = fill(ControlType.None, length(node_id)) + set_allocation_controlled!(control_type, node_id, parsed_parameters.control_mapping) + + Pump(; node_id, inflow_link = inflow_link.(Ref(graph), node_id), outflow_link = outflow_link.(Ref(graph), node_id), @@ -754,6 +757,7 @@ function Pump(db::DB, config::Config, graph::MetaGraph)::Pump parsed_parameters.min_upstream_level, parsed_parameters.max_downstream_level, parsed_parameters.control_mapping, + control_type, ) end @@ -801,6 +805,9 @@ function Outlet(db::DB, config::Config, graph::MetaGraph)::Outlet flow_rate_cache[Float64[], length(node_id)] .= [itp(0) for itp in parsed_parameters.flow_rate] + control_type = fill(ControlType.None, length(node_id)) + set_allocation_controlled!(control_type, node_id, parsed_parameters.control_mapping) + return Outlet(; node_id, inflow_link = inflow_link.(Ref(graph), node_id), @@ -813,6 +820,7 @@ function Outlet(db::DB, config::Config, graph::MetaGraph)::Outlet parsed_parameters.control_mapping, parsed_parameters.min_upstream_level, parsed_parameters.max_downstream_level, + control_type, ) end @@ -1878,7 +1886,7 @@ function Parameters(db::DB, config::Config)::Parameters ) collect_control_mappings!(p) - set_continuous_control_type!(p) + set_control_type!(p) set_listen_variable_refs!(p) set_discrete_controlled_variable_refs!(p) set_continuously_controlled_variable_refs!(p) diff --git a/core/src/solve.jl b/core/src/solve.jl index bcb23b8c7..08c35d3a2 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -36,7 +36,7 @@ function water_balance!(du::Vector, u::Vector, p::Parameters, t::Number)::Nothin t, current_low_storage_factor, current_level; - continuous_control_type = ContinuousControlType.Continuous, + control_type = ControlType.Continuous, ) # Compute PID control @@ -49,7 +49,7 @@ function water_balance!(du::Vector, u::Vector, p::Parameters, t::Number)::Nothin t, current_low_storage_factor, current_level; - continuous_control_type = ContinuousControlType.PID, + control_type = ControlType.PID, ) return nothing @@ -562,7 +562,7 @@ function formulate_flow!( t::Number, current_low_storage_factor::Vector, current_level::Vector, - continuous_control_type_::ContinuousControlType.T, + control_type_::ControlType.T, )::Nothing du_pump = view(du, p.state_ranges.pump) all_nodes_active = p.all_nodes_active[] @@ -577,7 +577,7 @@ function formulate_flow!( max_flow_rate, min_upstream_level, max_downstream_level, - continuous_control_type, + control_type, ) in zip( pump.node_id, pump.inflow_link, @@ -589,14 +589,13 @@ function formulate_flow!( pump.max_flow_rate, pump.min_upstream_level, pump.max_downstream_level, - pump.continuous_control_type, + pump.control_type, ) - if !(active || all_nodes_active) || - (continuous_control_type != continuous_control_type_) + if !(active || all_nodes_active) || (control_type != control_type_) continue end - flow_rate = if continuous_control_type == ContinuousControlType.None + flow_rate = if control_type == ControlType.None flow_rate_itp(t) else flow_rate_from_cache @@ -626,7 +625,7 @@ function formulate_flow!( t::Number, current_low_storage_factor::Vector, current_level::Vector, - continuous_control_type_::ContinuousControlType.T, + control_type_::ControlType.T, )::Nothing du_outlet = view(du, p.state_ranges.outlet) all_nodes_active = p.all_nodes_active[] @@ -639,7 +638,7 @@ function formulate_flow!( flow_rate_itp, min_flow_rate, max_flow_rate, - continuous_control_type, + control_type, min_upstream_level, max_downstream_level, ) in zip( @@ -651,16 +650,15 @@ function formulate_flow!( outlet.flow_rate, outlet.min_flow_rate, outlet.max_flow_rate, - outlet.continuous_control_type, + outlet.control_type, outlet.min_upstream_level, outlet.max_downstream_level, ) - if !(active || all_nodes_active) || - (continuous_control_type != continuous_control_type_) + if !(active || all_nodes_active) || (control_type != control_type_) continue end - flow_rate = if continuous_control_type == ContinuousControlType.None + flow_rate = if control_type == ControlType.None flow_rate_itp(t) else flow_rate_from_cache @@ -692,7 +690,7 @@ function formulate_flows!( t::Number, current_low_storage_factor::Vector, current_level::Vector; - continuous_control_type::ContinuousControlType.T = ContinuousControlType.None, + control_type::ControlType.T = ControlType.None, )::Nothing (; linear_resistance, @@ -703,15 +701,7 @@ function formulate_flows!( user_demand, ) = p - formulate_flow!( - du, - pump, - p, - t, - current_low_storage_factor, - current_level, - continuous_control_type, - ) + formulate_flow!(du, pump, p, t, current_low_storage_factor, current_level, control_type) formulate_flow!( du, outlet, @@ -719,10 +709,10 @@ function formulate_flows!( t, current_low_storage_factor, current_level, - continuous_control_type, + control_type, ) - if continuous_control_type == ContinuousControlType.None + if control_type == ControlType.None formulate_flow!( du, linear_resistance, diff --git a/core/src/util.jl b/core/src/util.jl index 7ab5fd345..19ada598a 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -400,20 +400,14 @@ function get_external_demand_priority_idx(p::Parameters, node_id::NodeID)::Int end """ -Set continuous_control_type for those pumps and outlets that are controlled by either +Set control_type for those pumps and outlets that are controlled by either PidControl or ContinuousControl """ -function set_continuous_control_type!(p::Parameters)::Nothing +function set_control_type!(p::Parameters)::Nothing (; continuous_control, pid_control) = p - errors = false - errors = set_continuous_control_type!( - p, - continuous_control.node_id, - ContinuousControlType.Continuous, - ) - errors |= - set_continuous_control_type!(p, pid_control.node_id, ContinuousControlType.PID) + errors = set_control_type!(p, continuous_control.node_id, ControlType.Continuous) + errors |= set_control_type!(p, pid_control.node_id, ControlType.PID) if errors error("Errors occurred when parsing ContinuousControl and PidControl connectivity") @@ -421,10 +415,10 @@ function set_continuous_control_type!(p::Parameters)::Nothing return nothing end -function set_continuous_control_type!( +function set_control_type!( p::Parameters, node_id::Vector{NodeID}, - continuous_control_type::ContinuousControlType.T, + control_type::ControlType.T, )::Bool (; graph, pump, outlet) = p errors = false @@ -432,9 +426,9 @@ function set_continuous_control_type!( for id in node_id id_controlled = only(outneighbor_labels_type(graph, id, LinkType.control)) if id_controlled.type == NodeType.Pump - pump.continuous_control_type[id_controlled.idx] = continuous_control_type + pump.control_type[id_controlled.idx] = control_type elseif id_controlled.type == NodeType.Outlet - outlet.continuous_control_type[id_controlled.idx] = continuous_control_type + outlet.control_type[id_controlled.idx] = control_type else errors = true @error "Only Pump and Outlet can be controlled by PidController, got $id_controlled" @@ -1146,3 +1140,20 @@ function get_timeseries_tstops( return tstops end + +""" +If for a certain node the "Ribasim.allocation" control state is defined, +make this the (initial) control type of that node +""" +function set_allocation_controlled!( + control_type::Vector{ControlType.T}, + node_id::Vector{NodeID}, + control_mapping::Dict{Tuple{NodeID, String}, ControlStateUpdate}, +)::Nothing + for id in node_id + if (id, "Ribasim.allocation") in keys(control_mapping) + control_type[id.idx] = ControlType.Allocation + end + end + return nothing +end diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index 673a72a93..6317d6455 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -6,6 +6,7 @@ import ribasim_testmodels from ribasim_testmodels.allocation import ( + allocation_control_model, allocation_example_model, allocation_training_model, bommelerwaard_model, @@ -71,6 +72,7 @@ from ribasim_testmodels.two_basin import two_basin_model __all__ = [ + "allocation_control_model", "allocation_training_model", "allocation_example_model", "backwater_model", diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index 9e3fa7e3d..ea259450d 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -859,7 +859,7 @@ def flow_demand_model() -> Model: return model -def linear_resistance_demand_model(): +def linear_resistance_demand_model() -> Model: """Small model with a FlowDemand for a node with a max flow rate.""" model = Model( starttime="2020-01-01 00:00:00", @@ -895,7 +895,7 @@ def linear_resistance_demand_model(): return model -def fair_distribution_model(): +def fair_distribution_model() -> Model: """See the behavior of allocation with few restrictions within the graph.""" model = Model( starttime="2020-01-01 00:00:00", @@ -994,7 +994,7 @@ def fair_distribution_model(): return model -def allocation_training_model(): +def allocation_training_model() -> Model: model = Model( starttime="2022-01-01", endtime="2023-01-01", @@ -1173,7 +1173,7 @@ def allocation_training_model(): return model -def bommelerwaard_model(): +def bommelerwaard_model() -> Model: model = Model( starttime="2016-01-01", endtime="2016-03-31", @@ -1368,7 +1368,7 @@ def get_node_data(node_id): return model -def cyclic_demand_model(): +def cyclic_demand_model() -> Model: """Create a model that has cyclic User- Flow- and LevelDemand.""" model = Model( starttime="2020-01-01", @@ -1442,3 +1442,50 @@ def cyclic_demand_model(): model.link.add(ud, bsn2) return model + + +def allocation_control_model(): + """Create a model that has a pump controlled by allocation.""" + model = Model( + starttime="2020-01-01", + endtime="2023-01-01", + crs="EPSG:28992", + allocation=Allocation(use_allocation=True), + ) + + lb = model.level_boundary.add( + Node(1, Point(0, 0), subnetwork_id=1), [level_boundary.Static(level=[1.0])] + ) + + out = model.outlet.add( + Node(2, Point(1, 0), subnetwork_id=1), + [outlet.Static(flow_rate=[0.0], control_state="Ribasim.allocation")], + ) + + bsn = model.basin.add( + Node(3, Point(2, 0), subnetwork_id=1), + [ + basin.State(level=[1.0]), + basin.Profile(level=[0.0, 1.0], area=[100.0, 100.0]), + ], + ) + + user = model.user_demand.add( + Node(4, Point(2, -1), subnetwork_id=1, cyclic_time=True), + [ + user_demand.Time( + time=["2020-01-01", "2020-06-01", "2021-01-01"], + demand=[0.0, 10.0, 0.0], + return_factor=0.0, + min_level=-5.0, + demand_priority=1, + ) + ], + ) + + model.link.add(lb, out) + model.link.add(out, bsn) + model.link.add(bsn, user) + model.link.add(user, bsn) + + return model