Code
using Ribasim
@@ -598,25 +602,25 @@
Min F[(Basin #2, UserDemand #3)]² + F[(Basin #5, UserDemand #6)]²
Subject to
- flow_conservation[TabulatedRatingCurve #7] : -F[(TabulatedRatingCurve #7, Terminal #8)] + F[(Basin #5, TabulatedRatingCurve #7)] = 0
- flow_conservation[Basin #2] : F[(UserDemand #3, Basin #2)] - F[(Basin #2, UserDemand #3)] - F[(Basin #2, LinearResistance #4)] + F[(LinearResistance #4, Basin #2)] + F[(FlowBoundary #1, Basin #2)] = 0
- flow_conservation[Basin #5] : F[(UserDemand #6, Basin #5)] + F[(LinearResistance #4, Basin #5)] - F[(Basin #5, LinearResistance #4)] - F[(Basin #5, UserDemand #6)] - F[(Basin #5, TabulatedRatingCurve #7)] = 0
- flow_conservation[LinearResistance #4] : F[(Basin #2, LinearResistance #4)] - F[(LinearResistance #4, Basin #2)] - F[(LinearResistance #4, Basin #5)] + F[(Basin #5, LinearResistance #4)] = 0
+ flow_conservation[Basin #5] : -F[(Basin #5, TabulatedRatingCurve #7)] + F[(LinearResistance #4, Basin #5)] - F[(Basin #5, LinearResistance #4)] + F[(UserDemand #6, Basin #5)] - F[(Basin #5, UserDemand #6)] = 0
flow_conservation[Terminal #8] : F[(TabulatedRatingCurve #7, Terminal #8)] = 0
+ flow_conservation[LinearResistance #4] : -F[(LinearResistance #4, Basin #5)] + F[(Basin #5, LinearResistance #4)] + F[(Basin #2, LinearResistance #4)] - F[(LinearResistance #4, Basin #2)] = 0
+ flow_conservation[Basin #2] : F[(FlowBoundary #1, Basin #2)] - F[(Basin #2, UserDemand #3)] + F[(UserDemand #3, Basin #2)] - F[(Basin #2, LinearResistance #4)] + F[(LinearResistance #4, Basin #2)] = 0
+ flow_conservation[TabulatedRatingCurve #7] : F[(Basin #5, TabulatedRatingCurve #7)] - F[(TabulatedRatingCurve #7, Terminal #8)] = 0
source[(FlowBoundary #1, Basin #2)] : F[(FlowBoundary #1, Basin #2)] ≤ 172800
source_user[UserDemand #6] : F[(UserDemand #6, Basin #5)] ≤ 0
source_user[UserDemand #3] : F[(UserDemand #3, Basin #2)] ≤ 0
+ F[(FlowBoundary #1, Basin #2)] ≥ 0
+ F[(Basin #5, TabulatedRatingCurve #7)] ≥ 0
+ F[(LinearResistance #4, Basin #5)] ≥ 0
+ F[(Basin #5, LinearResistance #4)] ≥ 0
F[(UserDemand #6, Basin #5)] ≥ 0
+ F[(Basin #2, UserDemand #3)] ≥ 0
F[(TabulatedRatingCurve #7, Terminal #8)] ≥ 0
+ F[(Basin #5, UserDemand #6)] ≥ 0
F[(UserDemand #3, Basin #2)] ≥ 0
- F[(Basin #2, UserDemand #3)] ≥ 0
F[(Basin #2, LinearResistance #4)] ≥ 0
F[(LinearResistance #4, Basin #2)] ≥ 0
- F[(FlowBoundary #1, Basin #2)] ≥ 0
- F[(LinearResistance #4, Basin #5)] ≥ 0
- F[(Basin #5, LinearResistance #4)] ≥ 0
- F[(Basin #5, UserDemand #6)] ≥ 0
- F[(Basin #5, TabulatedRatingCurve #7)] ≥ 0
- - -
- - diff --git a/guide/delwaq.html b/guide/delwaq.html index 2bb48005c..b11b2592a 100644 --- a/guide/delwaq.html +++ b/guide/delwaq.html @@ -134,7 +134,11 @@ Overview
- - + +Tutorials + +
- + How-to guides
-
@@ -181,12 +185,6 @@
- - -
- - diff --git a/guide/examples.html b/guide/examples.html index 3cd79ee12..ccf1c1bb2 100644 --- a/guide/examples.html +++ b/guide/examples.html @@ -65,7 +65,6 @@ - @@ -161,7 +160,11 @@ Overview
- - + +Tutorials + +
- + How-to guides
-
@@ -208,12 +211,6 @@
- - -
-
-
@@ -492,21 +489,14 @@
1 Basic model wit
1.1 Running a model
Now run the model. You can open a terminal and run it from there. For example, to run the basic model, input:
-From Python you can run it with:
--import subprocess -result = subprocess.run([cli_path, toml_path], capture_output=True, encoding="utf-8") -print(result.stderr) -result.check_returncode()
Where
-cli_path
is a string with either the full path to the Ribasim executable, liker"c:\ribasim_windows\ribasim"
, or just"ribasim"
in case you added theribasim_windows
folder to your PATH.The
print(result.stderr)
ensures you see the same logging and error messages that you would see in the terminal. Andresult.check_returncode()
will throw an error when the simulation was not successful.After running the model, read back the results:
-+df_basin = pd.read_feather(datadir / "basic/results/basin.arrow") -df_basin_wide = df_basin.pivot_table( - index="time", columns="node_id", values=["storage", "level"] -) -ax = df_basin_wide["level"].plot() -ax.set_ylabel("level [m]");
df_basin = pd.read_feather(datadir / "basic/results/basin.arrow") +df_basin_wide = df_basin.pivot_table( + index="time", columns="node_id", values=["storage", "level"] +) +ax = df_basin_wide["level"].plot() +ax.set_ylabel("level [m]");
-+df_flow = pd.read_feather(datadir / "basic/results/flow.arrow") -df_flow["edge"] = list(zip(df_flow.from_node_id, df_flow.to_node_id)) -df_flow["flow_m3d"] = df_flow.flow_rate * 86400 -ax = df_flow.pivot_table(index="time", columns="edge", values="flow_m3d").plot() -ax.legend(bbox_to_anchor=(1.3, 1), title="Edge") -ax.set_ylabel("flow [m³day⁻¹]");
df_flow = pd.read_feather(datadir / "basic/results/flow.arrow") +df_flow["edge"] = list(zip(df_flow.from_node_id, df_flow.to_node_id)) +df_flow["flow_m3d"] = df_flow.flow_rate * 86400 +ax = df_flow.pivot_table(index="time", columns="edge", values="flow_m3d").plot() +ax.legend(bbox_to_anchor=(1.3, 1), title="Edge") +ax.set_ylabel("flow [m³day⁻¹]");
@@ -537,42 +527,42 @@2 Model with disc
The model constructed below consists of a single basin which slowly drains trough a
TabulatedRatingCurve
, but is held within a range by two connected pumps. These two pumps together behave like a reversible pump. When pumping can be done in only one direction, and the other direction is only possible under gravity, use an Outlet for that direction.Setup the basins:
- +-+model.basin.add( - Node(1, Point(0.0, 0.0)), - [ - basin.Profile(area=[1000.0, 1000.0], level=[0.0, 1.0]), - basin.State(level=[20.0]), - basin.Time(time=["2020-01-01", "2020-07-01"], precipitation=[0.0, 3e-6]), - ], -)
model.basin.add( + Node(1, Point(0.0, 0.0)), + [ + basin.Profile(area=[1000.0, 1000.0], level=[0.0, 1.0]), + basin.State(level=[20.0]), + basin.Time(time=["2020-01-01", "2020-07-01"], precipitation=[0.0, 3e-6]), + ], +)
NodeData(node_id=1, node_type='Basin', geometry=<POINT (0 0)>)
Setup the discrete control:
-+model.discrete_control.add( - Node(7, Point(1.0, 0.0)), - [ - discrete_control.Variable( - compound_variable_id=1, - listen_node_id=1, - variable=["level"], - ), - discrete_control.Condition( - compound_variable_id=1, - # min, max - greater_than=[5.0, 15.0], - ), - discrete_control.Logic( - truth_state=["FF", "TF", "TT"], - control_state=["in", "none", "out"], - ), - ], -)
model.discrete_control.add( + Node(7, Point(1.0, 0.0)), + [ + discrete_control.Variable( + compound_variable_id=1, + listen_node_id=1, + variable=["level"], + ), + discrete_control.Condition( + compound_variable_id=1, + # min, max + greater_than=[5.0, 15.0], + ), + discrete_control.Logic( + truth_state=["FF", "TF", "TT"], + control_state=["in", "none", "out"], + ), + ], +)
NodeData(node_id=7, node_type='DiscreteControl', geometry=<POINT (1 0)>)
2 Model with disc
Setup the pump:
-+model.pump.add( - Node(2, Point(1.0, 1.0)), - [pump.Static(control_state=["none", "in", "out"], flow_rate=[0.0, 2e-3, 0.0])], -) -model.pump.add( - Node(3, Point(1.0, -1.0)), - [pump.Static(control_state=["none", "in", "out"], flow_rate=[0.0, 0.0, 2e-3])], -)
model.pump.add( + Node(2, Point(1.0, 1.0)), + [pump.Static(control_state=["none", "in", "out"], flow_rate=[0.0, 2e-3, 0.0])], +) +model.pump.add( + Node(3, Point(1.0, -1.0)), + [pump.Static(control_state=["none", "in", "out"], flow_rate=[0.0, 0.0, 2e-3])], +)
NodeData(node_id=3, node_type='Pump', geometry=<POINT (1 -1)>)
2 Model with disc
Setup the level boundary:
- +NodeData(node_id=4, node_type='LevelBoundary', geometry=<POINT (2 0)>)
Setup the rating curve:
-+model.tabulated_rating_curve.add( - Node(5, Point(-1.0, 0.0)), - [tabulated_rating_curve.Static(level=[2.0, 15.0], flow_rate=[0.0, 2e-3])], -)
model.tabulated_rating_curve.add( + Node(5, Point(-1.0, 0.0)), + [tabulated_rating_curve.Static(level=[2.0, 15.0], flow_rate=[0.0, 2e-3])], +)
NodeData(node_id=5, node_type='TabulatedRatingCurve', geometry=<POINT (-1 0)>)
Setup the terminal:
- +NodeData(node_id=6, node_type='Terminal', geometry=<POINT (-2 0)>)
Setup edges:
-+model.edge.add(model.basin[1], model.pump[3]) -model.edge.add(model.pump[3], model.level_boundary[4]) -model.edge.add(model.level_boundary[4], model.pump[2]) -model.edge.add(model.pump[2], model.basin[1]) -model.edge.add(model.basin[1], model.tabulated_rating_curve[5]) -model.edge.add(model.tabulated_rating_curve[5], model.terminal[6]) -model.edge.add(model.discrete_control[7], model.pump[2]) -model.edge.add(model.discrete_control[7], model.pump[3])
model.edge.add(model.basin[1], model.pump[3]) +model.edge.add(model.pump[3], model.level_boundary[4]) +model.edge.add(model.level_boundary[4], model.pump[2]) +model.edge.add(model.pump[2], model.basin[1]) +model.edge.add(model.basin[1], model.tabulated_rating_curve[5]) +model.edge.add(model.tabulated_rating_curve[5], model.terminal[6]) +model.edge.add(model.discrete_control[7], model.pump[2]) +model.edge.add(model.discrete_control[7], model.pump[3])
Let’s take a look at the model:
- +/home/runner/work/Ribasim/Ribasim/python/ribasim/ribasim/model.py:484: FutureWarning: @@ -680,35 +670,35 @@
2 Model with disc
Listen edges are plotted with a dashed line since they are not present in the “Edge / static” schema but only in the “Control / condition” schema.
- +PosixPath('data/level_range/ribasim.toml')
Now run the model (for running instructions see here). After running the model, read back the results:
-+df_basin = pd.read_feather(datadir / "level_range/results/basin.arrow") -df_basin_wide = df_basin.pivot_table( - index="time", columns="node_id", values=["storage", "level"] -) - -ax = df_basin_wide["level"].plot() - -greater_than = model.discrete_control.condition.df.greater_than - -ax.hlines( - greater_than, - df_basin.time[0], - df_basin.time.max(), - lw=1, - ls="--", - color="k", -) - -ax.set_yticks(greater_than, ["min", "max"]) -ax.set_ylabel("level") -plt.show()
df_basin = pd.read_feather(datadir / "level_range/results/basin.arrow") +df_basin_wide = df_basin.pivot_table( + index="time", columns="node_id", values=["storage", "level"] +) + +ax = df_basin_wide["level"].plot() + +greater_than = model.discrete_control.condition.df.greater_than + +ax.hlines( + greater_than, + df_basin.time[0], + df_basin.time.max(), + lw=1, + ls="--", + color="k", +) + +ax.set_yticks(greater_than, ["min", "max"]) +ax.set_ylabel("level") +plt.show()
@@ -723,95 +713,95 @@2 Model with disc
3 Model with PID control
Set up the model:
- +Setup the basins:
-+model.basin.add( - Node(2, Point(1.0, 0.0)), - [basin.Profile(area=[1000.0, 1000.0], level=[0.0, 1.0]), basin.State(level=[6.0])], -)
model.basin.add( + Node(2, Point(1.0, 0.0)), + [basin.Profile(area=[1000.0, 1000.0], level=[0.0, 1.0]), basin.State(level=[6.0])], +)
NodeData(node_id=2, node_type='Basin', geometry=<POINT (1 0)>)
Setup the pump:
-+model.pump.add( - Node(3, Point(2.0, 0.5)), - [pump.Static(flow_rate=[0.0])], # Will be overwritten by PID controller -)
model.pump.add( + Node(3, Point(2.0, 0.5)), + [pump.Static(flow_rate=[0.0])], # Will be overwritten by PID controller +)
NodeData(node_id=3, node_type='Pump', geometry=<POINT (2 0.5)>)
Setup the outlet:
-+model.outlet.add( - Node(6, Point(2.0, -0.5)), - [outlet.Static(flow_rate=[0.0])], # Will be overwritten by PID controller -)
model.outlet.add( + Node(6, Point(2.0, -0.5)), + [outlet.Static(flow_rate=[0.0])], # Will be overwritten by PID controller +)
NodeData(node_id=6, node_type='Outlet', geometry=<POINT (2 -0.5)>)
Setup flow boundary:
- +NodeData(node_id=1, node_type='FlowBoundary', geometry=<POINT (0 0)>)
Setup level boundary:
- +NodeData(node_id=4, node_type='LevelBoundary', geometry=<POINT (3 0)>)
Setup PID control:
-+for node, proportional, integral in [ - (Node(5, Point(1.5, 1.0)), -1e-3, -1e-7), - (Node(7, Point(1.5, -1.0)), 1e-3, 1e-7), -]: - pid_control_data = [ - pid_control.Time( - time=[ - "2020-01-01", - "2020-05-01", - "2020-07-01", - "2020-12-01", - ], - listen_node_id=2, - target=[5.0, 5.0, 7.5, 7.5], - proportional=proportional, - integral=integral, - derivative=0.0, - ) - ] - model.pid_control.add(node, pid_control_data)
for node, proportional, integral in [ + (Node(5, Point(1.5, 1.0)), -1e-3, -1e-7), + (Node(7, Point(1.5, -1.0)), 1e-3, 1e-7), +]: + pid_control_data = [ + pid_control.Time( + time=[ + "2020-01-01", + "2020-05-01", + "2020-07-01", + "2020-12-01", + ], + listen_node_id=2, + target=[5.0, 5.0, 7.5, 7.5], + proportional=proportional, + integral=integral, + derivative=0.0, + ) + ] + model.pid_control.add(node, pid_control_data)
Note that the coefficients for the pump and the outlet are equal in magnitude but opposite in sign. This way the pump and the outlet equally work towards the same goal, while having opposite effects on the controlled basin due to their connectivity to this basin.
Setup the edges:
-+model.edge.add(model.flow_boundary[1], model.basin[2]) -model.edge.add(model.basin[2], model.pump[3]) -model.edge.add(model.pump[3], model.level_boundary[4]) -model.edge.add(model.level_boundary[4], model.outlet[6]) -model.edge.add(model.outlet[6], model.basin[2]) -model.edge.add(model.pid_control[5], model.pump[3]) -model.edge.add(model.pid_control[7], model.outlet[6])
model.edge.add(model.flow_boundary[1], model.basin[2]) +model.edge.add(model.basin[2], model.pump[3]) +model.edge.add(model.pump[3], model.level_boundary[4]) +model.edge.add(model.level_boundary[4], model.outlet[6]) +model.edge.add(model.outlet[6], model.basin[2]) +model.edge.add(model.pid_control[5], model.pump[3]) +model.edge.add(model.pid_control[7], model.outlet[6])
Let’s take a look at the model:
- +/home/runner/work/Ribasim/Ribasim/python/ribasim/ribasim/model.py:469: FutureWarning: @@ -828,28 +818,28 @@
3 Model with PID
Write the model to a TOML and GeoPackage:
- +PosixPath('data/pid_control/ribasim.toml')
Now run the model (for running instructions see here). After running the model, read back the results:
-+from matplotlib.dates import date2num - -df_basin = pd.read_feather(datadir / "pid_control/results/basin.arrow") -df_basin_wide = df_basin.pivot_table( - index="time", columns="node_id", values=["storage", "level"] -) -ax = df_basin_wide["level"].plot() -ax.set_ylabel("level [m]") - -# Plot target level -level_demands = model.pid_control.time.df.target.to_numpy()[:4] -times = date2num(model.pid_control.time.df.time)[:4] -ax.plot(times, level_demands, color="k", ls=":", label="target level") -pass
from matplotlib.dates import date2num + +df_basin = pd.read_feather(datadir / "pid_control/results/basin.arrow") +df_basin_wide = df_basin.pivot_table( + index="time", columns="node_id", values=["storage", "level"] +) +ax = df_basin_wide["level"].plot() +ax.set_ylabel("level [m]") + +# Plot target level +level_demands = model.pid_control.time.df.target.to_numpy()[:4] +times = date2num(model.pid_control.time.df.time)[:4] +ax.plot(times, level_demands, color="k", ls=":", label="target level") +pass
@@ -863,108 +853,108 @@3 Model with PID
4 Model with allocation (user demand)
Setup a model:
- +Setup the basins:
-+basin_data = [ - basin.Profile(area=[300_000.0, 300_000.0], level=[0.0, 1.0]), - basin.State(level=[1.0]), -] - -model.basin.add( - Node(2, Point(1.0, 0.0), subnetwork_id=1), - basin_data, -) -model.basin.add( - Node(5, Point(3.0, 0.0), subnetwork_id=1), - basin_data, -)
basin_data = [ + basin.Profile(area=[300_000.0, 300_000.0], level=[0.0, 1.0]), + basin.State(level=[1.0]), +] + +model.basin.add( + Node(2, Point(1.0, 0.0), subnetwork_id=1), + basin_data, +) +model.basin.add( + Node(5, Point(3.0, 0.0), subnetwork_id=1), + basin_data, +)
NodeData(node_id=5, node_type='Basin', geometry=<POINT (3 0)>)
Setup the flow boundary:
-+model.flow_boundary.add( - Node(1, Point(0.0, 0.0), subnetwork_id=1), [flow_boundary.Static(flow_rate=[2.0])] -)
model.flow_boundary.add( + Node(1, Point(0.0, 0.0), subnetwork_id=1), [flow_boundary.Static(flow_rate=[2.0])] +)
NodeData(node_id=1, node_type='FlowBoundary', geometry=<POINT (0 0)>)
Setup the linear resistance:
-+model.linear_resistance.add( - Node(4, Point(2.0, 0.0), subnetwork_id=1), - [linear_resistance.Static(resistance=[0.06])], -)
model.linear_resistance.add( + Node(4, Point(2.0, 0.0), subnetwork_id=1), + [linear_resistance.Static(resistance=[0.06])], +)
NodeData(node_id=4, node_type='LinearResistance', geometry=<POINT (2 0)>)
Setup the tabulated rating curve:
-+model.tabulated_rating_curve.add( - Node(7, Point(4.0, 0.0), subnetwork_id=1), - [tabulated_rating_curve.Static(level=[0.0, 0.5, 1.0], flow_rate=[0.0, 0.0, 2.0])], -)
model.tabulated_rating_curve.add( + Node(7, Point(4.0, 0.0), subnetwork_id=1), + [tabulated_rating_curve.Static(level=[0.0, 0.5, 1.0], flow_rate=[0.0, 0.0, 2.0])], +)
NodeData(node_id=7, node_type='TabulatedRatingCurve', geometry=<POINT (4 0)>)
Setup the terminal:
- +NodeData(node_id=8, node_type='Terminal', geometry=<POINT (5 0)>)
Setup the users:
-+model.user_demand.add( - Node(6, Point(3.0, 1.0), subnetwork_id=1), - [ - user_demand.Static( - demand=[1.5], return_factor=[0.0], min_level=[-1.0], priority=[1] - ) - ], -) -model.user_demand.add( - Node(3, Point(1.0, 1.0), subnetwork_id=1), - [ - user_demand.Time( - demand=[0.0, 1.0, 1.2, 1.2], - return_factor=[0.0, 0.0, 0.0, 0.0], - min_level=[-1.0, -1.0, -1.0, -1.0], - priority=[1, 1, 2, 2], - time=2 * ["2020-01-01", "2020-01-20"], - ) - ], -)
model.user_demand.add( + Node(6, Point(3.0, 1.0), subnetwork_id=1), + [ + user_demand.Static( + demand=[1.5], return_factor=[0.0], min_level=[-1.0], priority=[1] + ) + ], +) +model.user_demand.add( + Node(3, Point(1.0, 1.0), subnetwork_id=1), + [ + user_demand.Time( + demand=[0.0, 1.0, 1.2, 1.2], + return_factor=[0.0, 0.0, 0.0, 0.0], + min_level=[-1.0, -1.0, -1.0, -1.0], + priority=[1, 1, 2, 2], + time=2 * ["2020-01-01", "2020-01-20"], + ) + ], +)
NodeData(node_id=3, node_type='UserDemand', geometry=<POINT (1 1)>)
Setup the allocation:
- +Setup the edges:
-+model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=1) -model.edge.add(model.basin[2], model.user_demand[3]) -model.edge.add(model.basin[2], model.linear_resistance[4]) -model.edge.add(model.linear_resistance[4], model.basin[5]) -model.edge.add(model.basin[5], model.user_demand[6]) -model.edge.add(model.basin[5], model.tabulated_rating_curve[7]) -model.edge.add(model.user_demand[3], model.basin[2]) -model.edge.add(model.user_demand[6], model.basin[5]) -model.edge.add(model.tabulated_rating_curve[7], model.terminal[8])
model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=1) +model.edge.add(model.basin[2], model.user_demand[3]) +model.edge.add(model.basin[2], model.linear_resistance[4]) +model.edge.add(model.linear_resistance[4], model.basin[5]) +model.edge.add(model.basin[5], model.user_demand[6]) +model.edge.add(model.basin[5], model.tabulated_rating_curve[7]) +model.edge.add(model.user_demand[3], model.basin[2]) +model.edge.add(model.user_demand[6], model.basin[5]) +model.edge.add(model.tabulated_rating_curve[7], model.terminal[8])
Let’s take a look at the model:
- +Write the model to a TOML and GeoPackage:
- +PosixPath('data/allocation_example/ribasim.toml')
Now run the model (for running instructions see here). After running the model, read back the results:
-+import matplotlib.ticker as plticker - -df_allocation = pd.read_feather(datadir / "allocation_example/results/allocation.arrow") -df_allocation_wide = df_allocation.pivot_table( - index="time", - columns=["node_type", "node_id", "priority"], - values=["demand", "allocated", "realized"], -) -df_allocation_wide = df_allocation_wide.loc[:, (df_allocation_wide != 0).any(axis=0)] - -fig, axs = plt.subplots(1, 3, figsize=(8, 5)) - -df_allocation_wide["demand"].plot(ax=axs[0], ls=":") -df_allocation_wide["allocated"].plot(ax=axs[1], ls="--") -df_allocation_wide.xs(1, level="priority", axis=1)["realized"].plot( - ax=axs[2], color=["C0", "C2", "C3"] -) - -fig.tight_layout() -loc = plticker.MultipleLocator(2) - -axs[0].set_ylabel("level [m]") - -for ax, title in zip(axs, ["Demand", "Allocated", "Abstracted"]): - ax.set_title(title) - ax.set_ylim(0.0, 1.6) - ax.xaxis.set_major_locator(loc)
import matplotlib.ticker as plticker + +df_allocation = pd.read_feather(datadir / "allocation_example/results/allocation.arrow") +df_allocation_wide = df_allocation.pivot_table( + index="time", + columns=["node_type", "node_id", "priority"], + values=["demand", "allocated", "realized"], +) +df_allocation_wide = df_allocation_wide.loc[:, (df_allocation_wide != 0).any(axis=0)] + +fig, axs = plt.subplots(1, 3, figsize=(8, 5)) + +df_allocation_wide["demand"].plot(ax=axs[0], ls=":") +df_allocation_wide["allocated"].plot(ax=axs[1], ls="--") +df_allocation_wide.xs(1, level="priority", axis=1)["realized"].plot( + ax=axs[2], color=["C0", "C2", "C3"] +) + +fig.tight_layout() +loc = plticker.MultipleLocator(2) + +axs[0].set_ylabel("level [m]") + +for ax, title in zip(axs, ["Demand", "Allocated", "Abstracted"]): + ax.set_title(title) + ax.set_ylim(0.0, 1.6) + ax.xaxis.set_major_locator(loc)
@@ -1025,14 +1015,14 @@4 Model with allo
Abstraction is accumulated over all priorities per user.
-+df_basin = pd.read_feather(datadir / "allocation_example/results/basin.arrow") -df_basin_wide = df_basin.pivot_table( - index="time", columns="node_id", values=["storage", "level"] -) - -ax = df_basin_wide["level"].plot() -ax.set_title("Basin levels") -ax.set_ylabel("level [m]")
df_basin = pd.read_feather(datadir / "allocation_example/results/basin.arrow") +df_basin_wide = df_basin.pivot_table( + index="time", columns="node_id", values=["storage", "level"] +) + +ax = df_basin_wide["level"].plot() +ax.set_title("Basin levels") +ax.set_ylabel("level [m]")
Text(0, 0.5, 'level [m]')
4 Model with allo
5 Model with allocation (basin supply/demand)
Setup a model:
- +Setup the basins:
-+basin_data = [ - basin.Profile(area=[1e3, 1e3], level=[0.0, 1.0]), - basin.State(level=[0.5]), -] -model.basin.add( - Node(2, Point(1.0, 0.0), subnetwork_id=2), - [ - *basin_data, - basin.Time( - time=["2020-01-01", "2020-01-16"], - drainage=[0.0, 0.0], - potential_evaporation=[0.0, 0.0], - infiltration=[0.0, 0.0], - precipitation=[1e-6, 0.0], - ), - ], -) -model.basin.add( - Node(5, Point(2.0, -1.0), subnetwork_id=2), - [ - *basin_data, - basin.Static( - drainage=[0.0], - potential_evaporation=[0.0], - infiltration=[0.0], - precipitation=[0.0], - ), - ], -)
basin_data = [ + basin.Profile(area=[1e3, 1e3], level=[0.0, 1.0]), + basin.State(level=[0.5]), +] +model.basin.add( + Node(2, Point(1.0, 0.0), subnetwork_id=2), + [ + *basin_data, + basin.Time( + time=["2020-01-01", "2020-01-16"], + drainage=[0.0, 0.0], + potential_evaporation=[0.0, 0.0], + infiltration=[0.0, 0.0], + precipitation=[1e-6, 0.0], + ), + ], +) +model.basin.add( + Node(5, Point(2.0, -1.0), subnetwork_id=2), + [ + *basin_data, + basin.Static( + drainage=[0.0], + potential_evaporation=[0.0], + infiltration=[0.0], + precipitation=[0.0], + ), + ], +)
NodeData(node_id=5, node_type='Basin', geometry=<POINT (2 -1)>)
Setup the flow boundary:
-+model.flow_boundary.add( - Node(1, Point(0.0, 0.0), subnetwork_id=2), [flow_boundary.Static(flow_rate=[1e-3])] -)
model.flow_boundary.add( + Node(1, Point(0.0, 0.0), subnetwork_id=2), [flow_boundary.Static(flow_rate=[1e-3])] +)
NodeData(node_id=1, node_type='FlowBoundary', geometry=<POINT (0 0)>)
Setup level demand:
-+model.level_demand.add( - Node(4, Point(1.0, -1.0), subnetwork_id=2), - [level_demand.Static(priority=[1], min_level=[1.0], max_level=[1.5])], -)
model.level_demand.add( + Node(4, Point(1.0, -1.0), subnetwork_id=2), + [level_demand.Static(priority=[1], min_level=[1.0], max_level=[1.5])], +)
NodeData(node_id=4, node_type='LevelDemand', geometry=<POINT (1 -1)>)
Setup the users:
-+model.user_demand.add( - Node(3, Point(2.0, 0.0), subnetwork_id=2), - [ - user_demand.Static( - priority=[2], demand=[1.5e-3], return_factor=[0.2], min_level=[0.2] - ) - ], -)
model.user_demand.add( + Node(3, Point(2.0, 0.0), subnetwork_id=2), + [ + user_demand.Static( + priority=[2], demand=[1.5e-3], return_factor=[0.2], min_level=[0.2] + ) + ], +)
NodeData(node_id=3, node_type='UserDemand', geometry=<POINT (2 0)>)
Setup the allocation:
- +Setup the edges:
- +Let’s take a look at the model:
- +Write the model to a TOML and GeoPackage:
- +PosixPath('data/level_demand/ribasim.toml')
Now run the model (for running instructions see here). After running the model, read back the results:
-+df_basin = pd.read_feather(datadir / "level_demand/results/basin.arrow") -df_basin = df_basin[df_basin.node_id == 2] -df_basin_wide = df_basin.pivot_table( - index="time", columns="node_id", values=["storage", "level"] -) -ax = df_basin_wide["level"].plot(ylabel="level [m]")
df_basin = pd.read_feather(datadir / "level_demand/results/basin.arrow") +df_basin = df_basin[df_basin.node_id == 2] +df_basin_wide = df_basin.pivot_table( + index="time", columns="node_id", values=["storage", "level"] +) +ax = df_basin_wide["level"].plot(ylabel="level [m]")
@@ -1187,98 +1177,98 @@6 Guidance of mod
The polder system is composed of a sequence of level controlled polder basins with weirs inbetween each basin and an inlet and outlet to main system
- +All the polder basins are exposed to time varying forcings (precipitation, evaporation, drainage, infiltration) to mimic situations of water excess and water shortage.
In case of water excess, a pump in the most downstream polder will need to pump the surplus water to the main water system. In case of water shortage, an inlet at the most upstream polder will need to bring water into the cascase of polders. The main water system acts as a water source.
Model approach: All polder basins as well as the main water system are modelled with basin nodes. To let the system experience all 4 excess/shortage situation, forcing time series are made in a way that is adapting to them. Overall, assume that in one year, the system will experience precipitation (situation 1) in winter and early spring, precipitation shortage (situation 2) from late spring until early autumn. During situation 2, polder basin 4 will experience additional seepage (compoensating its shortage), and later polder basin 3 will also receive more seepage.
Setting up the basins:
-+time = pd.date_range(model.starttime, model.endtime) -day_of_year = time.day_of_year.to_numpy() - -precipitation = np.zeros(day_of_year.size) -precipitation[0:90] = 1.72e-8 -precipitation[330:366] = 1.72e-8 - -evaporation = np.zeros(day_of_year.size) -evaporation[130:270] = 2.87e-8 - -drainage = np.zeros(day_of_year.size) -drainage[120:270] = 0.4 * 2.87e-8 -drainage_3 = drainage.copy() -drainage_3[210:240] = 17 * 2.87e-8 -drainage_4 = drainage.copy() -drainage_4[160:240] = 13 * 2.87e-8 - -infiltration = np.zeros(day_of_year.size) -infiltration[0:90] = 5e-8 - -polder_profile = basin.Profile(area=[100, 100], level=[0.0, 3.0]) - -basin_time = [ - basin.Time( - time=pd.date_range(model.starttime, model.endtime), - drainage=drainage, - potential_evaporation=evaporation, - infiltration=0.0, - precipitation=precipitation, - ), -] - -basin_time4 = [ - basin.Time( - time=pd.date_range(model.starttime, model.endtime), - drainage=drainage_4, - potential_evaporation=evaporation, - infiltration=0.0, - precipitation=precipitation, - ), -] -basin_time3 = [ - basin.Time( - time=pd.date_range(model.starttime, model.endtime), - drainage=drainage_3, - potential_evaporation=evaporation, - infiltration=0.0, - precipitation=precipitation, - ), -] - -model.basin.add( - Node(1, Point(2.0, 0.0)), - [ - basin.State(level=[2.5]), - basin.Profile(area=[1000, 1000], level=[0.0, 3.0]), - basin.Time( - time=pd.date_range(model.starttime, model.endtime), - drainage=0.0, - potential_evaporation=0.0, - infiltration=0.0, - precipitation=0.0, - ), - ], -) -model.basin.add( - Node(4, Point(0.0, -2.0)), - [basin.State(level=[1.5]), polder_profile, *basin_time], -) -model.basin.add( - Node(6, Point(0.0, -4.0)), - [basin.State(level=[1.0]), polder_profile, *basin_time], -) -model.basin.add( - Node(8, Point(2.0, -4.0)), - [basin.State(level=[1.5]), polder_profile, *basin_time3], -) -model.basin.add( - Node(10, Point(4.0, -4.0)), - [basin.State(level=[1.3]), polder_profile, *basin_time4], -) -model.basin.add( - Node(12, Point(4.0, -2.0)), - [basin.State(level=[0.1]), polder_profile, *basin_time], -)
time = pd.date_range(model.starttime, model.endtime) +day_of_year = time.day_of_year.to_numpy() + +precipitation = np.zeros(day_of_year.size) +precipitation[0:90] = 1.72e-8 +precipitation[330:366] = 1.72e-8 + +evaporation = np.zeros(day_of_year.size) +evaporation[130:270] = 2.87e-8 + +drainage = np.zeros(day_of_year.size) +drainage[120:270] = 0.4 * 2.87e-8 +drainage_3 = drainage.copy() +drainage_3[210:240] = 17 * 2.87e-8 +drainage_4 = drainage.copy() +drainage_4[160:240] = 13 * 2.87e-8 + +infiltration = np.zeros(day_of_year.size) +infiltration[0:90] = 5e-8 + +polder_profile = basin.Profile(area=[100, 100], level=[0.0, 3.0]) + +basin_time = [ + basin.Time( + time=pd.date_range(model.starttime, model.endtime), + drainage=drainage, + potential_evaporation=evaporation, + infiltration=0.0, + precipitation=precipitation, + ), +] + +basin_time4 = [ + basin.Time( + time=pd.date_range(model.starttime, model.endtime), + drainage=drainage_4, + potential_evaporation=evaporation, + infiltration=0.0, + precipitation=precipitation, + ), +] +basin_time3 = [ + basin.Time( + time=pd.date_range(model.starttime, model.endtime), + drainage=drainage_3, + potential_evaporation=evaporation, + infiltration=0.0, + precipitation=precipitation, + ), +] + +model.basin.add( + Node(1, Point(2.0, 0.0)), + [ + basin.State(level=[2.5]), + basin.Profile(area=[1000, 1000], level=[0.0, 3.0]), + basin.Time( + time=pd.date_range(model.starttime, model.endtime), + drainage=0.0, + potential_evaporation=0.0, + infiltration=0.0, + precipitation=0.0, + ), + ], +) +model.basin.add( + Node(4, Point(0.0, -2.0)), + [basin.State(level=[1.5]), polder_profile, *basin_time], +) +model.basin.add( + Node(6, Point(0.0, -4.0)), + [basin.State(level=[1.0]), polder_profile, *basin_time], +) +model.basin.add( + Node(8, Point(2.0, -4.0)), + [basin.State(level=[1.5]), polder_profile, *basin_time3], +) +model.basin.add( + Node(10, Point(4.0, -4.0)), + [basin.State(level=[1.3]), polder_profile, *basin_time4], +) +model.basin.add( + Node(12, Point(4.0, -2.0)), + [basin.State(level=[0.1]), polder_profile, *basin_time], +)
NodeData(node_id=12, node_type='Basin', geometry=<POINT (4 -2)>)
6 Guidance of mod
Since the water level of the main system is at 2.5 meter above the reference level a pump is needed to remove the water from polder basin 5.
Setup the pumps:
- +NodeData(node_id=13, node_type='Pump', geometry=<POINT (4 -1)>)
6 Guidance of mod
According to the description of situation 1 and 2, the water in one polder basin needs to be able to flow to the downstream basin if the current basin has too much water (i.e. the water level is above the setpoint) or if the downstream basin is below setpoint and needs more water. This could be modelled with an uncontrolled TabulatedRatingCurve node with Q=0 at the setpoint level (and Q rising when the level rises above setpoint) , or with an Outlet node where the
min_upstream_level
is specified at or just below the setpoint. In this example, we’ve chosen for the Outlet where we specify the minimum upstream level 5 cm below the setpoint. For example: the Outlet of polder basin 1 (node 4) is specified with a minimum upstream level of 1.95 meter.Setup the outlets:
-+# Set up outlet -model.outlet.add( - Node(2, Point(0.0, -1.0)), - [outlet.Static(flow_rate=[2 * 0.5 / 3600], min_upstream_level=[0.0])], -) -model.outlet.add( - Node(5, Point(0.0, -3.0)), - [outlet.Static(flow_rate=[0.5 / 3600], min_upstream_level=[1.95])], -) -model.outlet.add( - Node(7, Point(1.0, -4.0)), - [outlet.Static(flow_rate=[0.5 / 3600], min_upstream_level=[1.45])], -) -model.outlet.add( - Node(9, Point(3.0, -4.0)), - [outlet.Static(flow_rate=[0.5 / 3600], min_upstream_level=[0.95])], -) -model.outlet.add( - Node(11, Point(4.0, -3.0)), - [outlet.Static(flow_rate=[0.5 / 3600], min_upstream_level=[0.45])], -)
# Set up outlet +model.outlet.add( + Node(2, Point(0.0, -1.0)), + [outlet.Static(flow_rate=[2 * 0.5 / 3600], min_upstream_level=[0.0])], +) +model.outlet.add( + Node(5, Point(0.0, -3.0)), + [outlet.Static(flow_rate=[0.5 / 3600], min_upstream_level=[1.95])], +) +model.outlet.add( + Node(7, Point(1.0, -4.0)), + [outlet.Static(flow_rate=[0.5 / 3600], min_upstream_level=[1.45])], +) +model.outlet.add( + Node(9, Point(3.0, -4.0)), + [outlet.Static(flow_rate=[0.5 / 3600], min_upstream_level=[0.95])], +) +model.outlet.add( + Node(11, Point(4.0, -3.0)), + [outlet.Static(flow_rate=[0.5 / 3600], min_upstream_level=[0.45])], +)
NodeData(node_id=11, node_type='Outlet', geometry=<POINT (4 -3)>)
6 Guidance of mod
When using Outlets as connecting nodes, the flow over the Outlet needs to be controlled to maintain the water level at the setpoint. For this purpose we introduce local PidControllers, where the targets of the PidControllers are set to the setpoints. Disadvantage of this local control approach is the delay that is introduced to transport the ‘basin X has a shortage’ message upstream through the cascade to the inlet. Current functionality does not offer the capability for PidControl to take multiple observations into account when controlling the inlet. Combining multiple observations in one control is feasible with DiscreteControl. This could be an alternative approach to controlling the inlet for the cascading water system.
Setup the PID control:
-+pid_control_data = { - "proportional": [0.05], - "integral": [0.00], - "derivative": [0.0], -} -model.pid_control.add( - Node(3, Point(-1.0, -1.0)), - [pid_control.Static(listen_node_id=[4], target=[2.0], **pid_control_data)], -) -model.pid_control.add( - Node(14, Point(-1.0, -3.0)), - [pid_control.Static(listen_node_id=[6], target=[1.5], **pid_control_data)], -) -model.pid_control.add( - Node(15, Point(1.0, -3.0)), - [pid_control.Static(listen_node_id=[8], target=[1.0], **pid_control_data)], -) -model.pid_control.add( - Node(16, Point(3.0, -3.0)), - [pid_control.Static(listen_node_id=[10], target=[0.5], **pid_control_data)], -)
pid_control_data = { + "proportional": [0.05], + "integral": [0.00], + "derivative": [0.0], +} +model.pid_control.add( + Node(3, Point(-1.0, -1.0)), + [pid_control.Static(listen_node_id=[4], target=[2.0], **pid_control_data)], +) +model.pid_control.add( + Node(14, Point(-1.0, -3.0)), + [pid_control.Static(listen_node_id=[6], target=[1.5], **pid_control_data)], +) +model.pid_control.add( + Node(15, Point(1.0, -3.0)), + [pid_control.Static(listen_node_id=[8], target=[1.0], **pid_control_data)], +) +model.pid_control.add( + Node(16, Point(3.0, -3.0)), + [pid_control.Static(listen_node_id=[10], target=[0.5], **pid_control_data)], +)
NodeData(node_id=16, node_type='PidControl', geometry=<POINT (3 -3)>)
Setup the edges:
-+model.edge.add(model.basin[1], model.outlet[2]) -model.edge.add(model.pid_control[3], model.outlet[2]) -model.edge.add(model.outlet[2], model.basin[4]) -model.edge.add(model.basin[4], model.outlet[5]) -model.edge.add(model.outlet[5], model.basin[6]) -model.edge.add(model.basin[6], model.outlet[7]) -model.edge.add(model.outlet[7], model.basin[8]) -model.edge.add(model.basin[8], model.outlet[9]) -model.edge.add(model.outlet[9], model.basin[10]) -model.edge.add(model.basin[10], model.outlet[11]) -model.edge.add(model.outlet[11], model.basin[12]) -model.edge.add(model.basin[12], model.pump[13]) -model.edge.add(model.pump[13], model.basin[1]) -model.edge.add(model.pid_control[14], model.outlet[5]) -model.edge.add(model.pid_control[15], model.outlet[7]) -model.edge.add(model.pid_control[16], model.outlet[9])
model.edge.add(model.basin[1], model.outlet[2]) +model.edge.add(model.pid_control[3], model.outlet[2]) +model.edge.add(model.outlet[2], model.basin[4]) +model.edge.add(model.basin[4], model.outlet[5]) +model.edge.add(model.outlet[5], model.basin[6]) +model.edge.add(model.basin[6], model.outlet[7]) +model.edge.add(model.outlet[7], model.basin[8]) +model.edge.add(model.basin[8], model.outlet[9]) +model.edge.add(model.outlet[9], model.basin[10]) +model.edge.add(model.basin[10], model.outlet[11]) +model.edge.add(model.outlet[11], model.basin[12]) +model.edge.add(model.basin[12], model.pump[13]) +model.edge.add(model.pump[13], model.basin[1]) +model.edge.add(model.pid_control[14], model.outlet[5]) +model.edge.add(model.pid_control[15], model.outlet[7]) +model.edge.add(model.pid_control[16], model.outlet[9])
To plot the model
- +/home/runner/work/Ribasim/Ribasim/python/ribasim/ribasim/model.py:469: FutureWarning: @@ -1389,29 +1379,29 @@
6 Guidance of mod
Write the model to a TOML file and run it in the Julia.
- +PosixPath('data/local_pidcontrolled_cascade/ribasim.toml')
Now run the model (for running instructions see here). After running the model, read back the result to plot the flow of each polder basin.
-+datadir_flow = datadir / "local_pidcontrolled_cascade/results/flow.arrow" -df_flow = pd.read_feather(datadir_flow) -df_flow["edge"] = list(zip(df_flow.from_node_id, df_flow.to_node_id)) -df_flow["flow_m3d"] = df_flow.flow_rate * 86400 - -df_pivot = df_flow.pivot_table(index="time", columns="edge", values="flow_m3d")
datadir_flow = datadir / "local_pidcontrolled_cascade/results/flow.arrow" +df_flow = pd.read_feather(datadir_flow) +df_flow["edge"] = list(zip(df_flow.from_node_id, df_flow.to_node_id)) +df_flow["flow_m3d"] = df_flow.flow_rate * 86400 + +df_pivot = df_flow.pivot_table(index="time", columns="edge", values="flow_m3d")
Below graphs show the flow exchanged with the mainsystem (i.e. the inlet and the pump), and the flow of weirs inbetween the polder basins.
-+df_input = df_pivot.loc[:, [(1, 2), (13, 1)]] -ax = df_input.plot(ylim=[-1.0, 20.0]) -ax.set_ylabel("flow [m³day⁻¹]") -df_weirs = df_pivot.loc[:, [(4, 5), (6, 7), (8, 9), (10, 11)]] -ax = df_weirs.plot(ylim=[-1.0, 15.0]) -ax.set_ylabel("flow [m³day⁻¹]");
df_input = df_pivot.loc[:, [(1, 2), (13, 1)]] +ax = df_input.plot(ylim=[-1.0, 20.0]) +ax.set_ylabel("flow [m³day⁻¹]") +df_weirs = df_pivot.loc[:, [(4, 5), (6, 7), (8, 9), (10, 11)]] +ax = df_weirs.plot(ylim=[-1.0, 15.0]) +ax.set_ylabel("flow [m³day⁻¹]");
Below graph shows the vertical flux on each basin.
-+datadir_basin = datadir / "local_pidcontrolled_cascade/results/basin.arrow" -df_basin = pd.read_feather(datadir_basin) -df_basin["vertical_flux"] = ( - df_basin["precipitation"] - - df_basin["evaporation"] - + df_basin["drainage"] - + df_basin["infiltration"] -) -df_basin_wide = df_basin.pivot_table( - index="time", columns="node_id", values=["storage", "level", "vertical_flux"] -) -df_basin_wide["vertical_flux"] *= 86400 -ax = df_basin_wide["vertical_flux"].plot() -ax.set_ylabel("vertical flux [m³day⁻¹]");
datadir_basin = datadir / "local_pidcontrolled_cascade/results/basin.arrow" +df_basin = pd.read_feather(datadir_basin) +df_basin["vertical_flux"] = ( + df_basin["precipitation"] + - df_basin["evaporation"] + + df_basin["drainage"] + + df_basin["infiltration"] +) +df_basin_wide = df_basin.pivot_table( + index="time", columns="node_id", values=["storage", "level", "vertical_flux"] +) +df_basin_wide["vertical_flux"] *= 86400 +ax = df_basin_wide["vertical_flux"].plot() +ax.set_ylabel("vertical flux [m³day⁻¹]");
In the following graph, the water level of basins are shown. The five polder basins are given starting levels that are different from their setpoints. It can be observed that in the beginning, the water level are changing and approaching to the set points. Later when the water levels are stable, they will not be affected by the forcing.
- +@@ -1467,91 +1457,91 @@6 Guidance of mod
7 Model with continuous control
- +Set up the transient level boundary:
-+model.level_boundary.add( - Node(1, Point(0, 0)), - [ - level_boundary.Time( - time=pd.date_range(start="2020-01-01", end="2021-01-01", periods=100), - level=6.0 + np.sin(np.linspace(0, 6 * np.pi, 100)), - ) - ], -)
model.level_boundary.add( + Node(1, Point(0, 0)), + [ + level_boundary.Time( + time=pd.date_range(start="2020-01-01", end="2021-01-01", periods=100), + level=6.0 + np.sin(np.linspace(0, 6 * np.pi, 100)), + ) + ], +)
NodeData(node_id=1, node_type='LevelBoundary', geometry=<POINT (0 0)>)
Set up the linear resistance:
-+model.linear_resistance.add( - Node(2, Point(1, 0)), [linear_resistance.Static(resistance=[10.0])] -)
model.linear_resistance.add( + Node(2, Point(1, 0)), [linear_resistance.Static(resistance=[10.0])] +)
NodeData(node_id=2, node_type='LinearResistance', geometry=<POINT (1 0)>)
Set up the basin:
-+model.basin.add( - Node(3, Point(2, 0)), - [ - basin.Profile(area=10000.0, level=[0.0, 1.0]), - basin.State(level=[10.0]), - ], -)
model.basin.add( + Node(3, Point(2, 0)), + [ + basin.Profile(area=10000.0, level=[0.0, 1.0]), + basin.State(level=[10.0]), + ], +)
NodeData(node_id=3, node_type='Basin', geometry=<POINT (2 0)>)
Set up the outlets:
-+model.outlet.add(Node(4, Point(3, 1)), [outlet.Static(flow_rate=[1.0])]) -model.outlet.add(Node(5, Point(3, -1)), [outlet.Static(flow_rate=[1.0])])
model.outlet.add(Node(4, Point(3, 1)), [outlet.Static(flow_rate=[1.0])]) +model.outlet.add(Node(5, Point(3, -1)), [outlet.Static(flow_rate=[1.0])])
NodeData(node_id=5, node_type='Outlet', geometry=<POINT (3 -1)>)
Set up the terminals:
- +NodeData(node_id=7, node_type='Terminal', geometry=<POINT (4 -1)>)
Set up the continuous control:
-+model.continuous_control.add( - Node(8, Point(2, 1)), - [ - continuous_control.Variable( - listen_node_id=[2], - variable="flow_rate", - ), - continuous_control.Function( - input=[0.0, 1.0], - output=[0.0, 0.6], - controlled_variable="flow_rate", - ), - ], -) -model.continuous_control.add( - Node(9, Point(2, -1)), - [ - continuous_control.Variable( - listen_node_id=[2], - variable="flow_rate", - ), - continuous_control.Function( - input=[0.0, 1.0], - output=[0.0, 0.4], - controlled_variable="flow_rate", - ), - ], -)
model.continuous_control.add( + Node(8, Point(2, 1)), + [ + continuous_control.Variable( + listen_node_id=[2], + variable="flow_rate", + ), + continuous_control.Function( + input=[0.0, 1.0], + output=[0.0, 0.6], + controlled_variable="flow_rate", + ), + ], +) +model.continuous_control.add( + Node(9, Point(2, -1)), + [ + continuous_control.Variable( + listen_node_id=[2], + variable="flow_rate", + ), + continuous_control.Function( + input=[0.0, 1.0], + output=[0.0, 0.4], + controlled_variable="flow_rate", + ), + ], +)
NodeData(node_id=9, node_type='ContinuousControl', geometry=<POINT (2 -1)>)
7 Model with cont
- A
ContinuousControl
node with ID 2, which listens to the flow rate of theLinearResistance
node with ID 1, puts that through the function \(f(x) = \max(0, 0.4x)\), and assigns the result to the flow rate of the node thisContinuousControl
node is controlling, which is defined by a (control) edge.-+model.edge.add(model.level_boundary[1], model.linear_resistance[2]) -model.edge.add(model.linear_resistance[2], model.basin[3]) -model.edge.add(model.basin[3], model.outlet[4]) -model.edge.add(model.basin[3], model.outlet[5]) -model.edge.add(model.outlet[4], model.terminal[6]) -model.edge.add(model.outlet[5], model.terminal[7]) - -# Define which node is controlled by each continuous control node -model.edge.add(model.continuous_control[8], model.outlet[4]) -model.edge.add(model.continuous_control[9], model.outlet[5])
model.edge.add(model.level_boundary[1], model.linear_resistance[2]) +model.edge.add(model.linear_resistance[2], model.basin[3]) +model.edge.add(model.basin[3], model.outlet[4]) +model.edge.add(model.basin[3], model.outlet[5]) +model.edge.add(model.outlet[4], model.terminal[6]) +model.edge.add(model.outlet[5], model.terminal[7]) + +# Define which node is controlled by each continuous control node +model.edge.add(model.continuous_control[8], model.outlet[4]) +model.edge.add(model.continuous_control[9], model.outlet[5])
Let’s take a look at the model:
- +/home/runner/work/Ribasim/Ribasim/python/ribasim/ribasim/model.py:484: FutureWarning: @@ -1593,35 +1583,35 @@
7 Model with cont
With this setup we want to split the flow coming into the basin into a 60% - 40% ratio.
Write the model to a TOML and GeoPackage:
- +PosixPath('data/outlet_continuous_control/ribasim.toml')
Now run the model (for running instructions see here). After running the model, read back the results:
-+df_flow = pd.read_feather(datadir / "outlet_continuous_control/results/flow.arrow") -fig, ax = plt.subplots() - - -def plot_edge_flow(from_node_type, from_node_id, to_node_type, to_node_id): - df_flow_filtered = df_flow[ - (df_flow["from_node_id"] == from_node_id) - & (df_flow["to_node_id"] == to_node_id) - ] - df_flow_filtered.plot( - x="time", - y="flow_rate", - ax=ax, - label=f"{from_node_type} #{from_node_id} → {to_node_type} #{to_node_id}", - ) - - -plot_edge_flow("LinearResistance", 1, "Basin", 1) -plot_edge_flow("Basin", 1, "Outlet", 1) -plot_edge_flow("Basin", 1, "Outlet", 2) -ax.set_ylabel("flow [m³s⁻¹]");
df_flow = pd.read_feather(datadir / "outlet_continuous_control/results/flow.arrow") +fig, ax = plt.subplots() + + +def plot_edge_flow(from_node_type, from_node_id, to_node_type, to_node_id): + df_flow_filtered = df_flow[ + (df_flow["from_node_id"] == from_node_id) + & (df_flow["to_node_id"] == to_node_id) + ] + df_flow_filtered.plot( + x="time", + y="flow_rate", + ax=ax, + label=f"{from_node_type} #{from_node_id} → {to_node_type} #{to_node_id}", + ) + + +plot_edge_flow("LinearResistance", 1, "Basin", 1) +plot_edge_flow("Basin", 1, "Outlet", 1) +plot_edge_flow("Basin", 1, "Outlet", 2) +ax.set_ylabel("flow [m³s⁻¹]");
@@ -2057,9 +2047,6 @@7 Model with cont