diff --git a/core/src/solve.jl b/core/src/solve.jl index a6b3cd791..b3fa3c34b 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -154,14 +154,16 @@ function continuous_control!( factor_outlet = 1.0 - # No flow of outlet if source level is lower than target level + # No flow out of outlet if source level is lower than reference level if !(src_level === nothing || dst_level === nothing) Δlevel = src_level - dst_level factor_outlet *= reduction_factor(Δlevel, 0.1) end - # No flow out outlet if source level is lower than minimum crest level + # No flow out of outlet if source level is lower than minimum crest level if src_level !== nothing + controlled_node_idx = findsorted(outlet.node_id, controlled_node_id) + factor_outlet *= reduction_factor( src_level - outlet.min_crest_level[controlled_node_idx], 0.1, diff --git a/docs/python/examples.ipynb b/docs/python/examples.ipynb index 27f3e3727..315ce4777 100644 --- a/docs/python/examples.ipynb +++ b/docs/python/examples.ipynb @@ -1677,7 +1677,436 @@ { "cell_type": "markdown", "metadata": {}, - "source": [] + "source": [ + "# Guidance of modelling a cascade of polder basins" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Situation description**: This example shows how to make a model for a given practical water system, which consists of a cascade of level control polder basins with inlet and outlet to the main systems. Note that alternative model layouts are feasible for the same water system, each having its positive items and drawbacks." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "![](image/cascade-polder.jpg)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "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" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model = Model(starttime=\"2020-01-01\", endtime=\"2021-01-01\", crs=\"EPSG:28992\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "All the polder basins are exposed to time varying forcings (precipitation, evaporation, drainage, infiltration) to mimic situations of water excess and water shortage.\n", + "\n", + "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.\n", + "\n", + "**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." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setting up the basins:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "time = pd.date_range(model.starttime, model.endtime)\n", + "day_of_year = time.day_of_year.to_numpy()\n", + "\n", + "precipitation = np.zeros(day_of_year.size)\n", + "precipitation[0:90] = 1.72e-8\n", + "precipitation[330:366] = 1.72e-8\n", + "\n", + "evaporation = np.zeros(day_of_year.size)\n", + "evaporation[130:270] = 2.87e-8\n", + "\n", + "drainage = np.zeros(day_of_year.size)\n", + "drainage[120:270] = 0.4 * 2.87e-8\n", + "drainage_3 = drainage.copy()\n", + "drainage_3[210:240] = 17 * 2.87e-8\n", + "drainage_4 = drainage.copy()\n", + "drainage_4[160:240] = 13 * 2.87e-8\n", + "\n", + "infiltration = np.zeros(day_of_year.size)\n", + "infiltration[0:90] = 5e-8\n", + "\n", + "polder_profile = basin.Profile(area=[100, 100], level=[0.0, 3.0])\n", + "\n", + "basin_time = [\n", + " basin.Time(\n", + " time=pd.date_range(model.starttime, model.endtime),\n", + " drainage=drainage,\n", + " potential_evaporation=evaporation,\n", + " infiltration=0.0,\n", + " precipitation=precipitation,\n", + " urban_runoff=0.0,\n", + " ),\n", + "]\n", + "\n", + "basin_time4 = [\n", + " basin.Time(\n", + " time=pd.date_range(model.starttime, model.endtime),\n", + " drainage=drainage_4,\n", + " potential_evaporation=evaporation,\n", + " infiltration=0.0,\n", + " precipitation=precipitation,\n", + " urban_runoff=0.0,\n", + " ),\n", + "]\n", + "basin_time3 = [\n", + " basin.Time(\n", + " time=pd.date_range(model.starttime, model.endtime),\n", + " drainage=drainage_3,\n", + " potential_evaporation=evaporation,\n", + " infiltration=0.0,\n", + " precipitation=precipitation,\n", + " urban_runoff=0.0,\n", + " ),\n", + "]\n", + "\n", + "model.basin.add(\n", + " Node(1, Point(2.0, 0.0)),\n", + " [\n", + " basin.State(level=[2.5]),\n", + " basin.Profile(area=[1000, 1000], level=[0.0, 3.0]),\n", + " basin.Time(\n", + " time=pd.date_range(model.starttime, model.endtime),\n", + " drainage=0.0,\n", + " potential_evaporation=0.0,\n", + " infiltration=0.0,\n", + " precipitation=0.0,\n", + " urban_runoff=0.0,\n", + " ),\n", + " ],\n", + ")\n", + "model.basin.add(\n", + " Node(4, Point(0.0, -2.0)),\n", + " [basin.State(level=[1.5]), polder_profile, *basin_time],\n", + ")\n", + "model.basin.add(\n", + " Node(6, Point(0.0, -4.0)),\n", + " [basin.State(level=[1.0]), polder_profile, *basin_time],\n", + ")\n", + "model.basin.add(\n", + " Node(8, Point(2.0, -4.0)),\n", + " [basin.State(level=[1.5]), polder_profile, *basin_time3],\n", + ")\n", + "model.basin.add(\n", + " Node(10, Point(4.0, -4.0)),\n", + " [basin.State(level=[1.3]), polder_profile, *basin_time4],\n", + ")\n", + "model.basin.add(\n", + " Node(12, Point(4.0, -2.0)),\n", + " [basin.State(level=[0.1]), polder_profile, *basin_time],\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After all the basins are defined the connecting component inbetween the basins needs to be determined. For polder basin 5 (node 12), the water level needs to be maintain at 0.0 meter. This means that either there should be no water in this basin, or the basin bottom is lower than the reference level, and the water level should be maintained at the reference level.\n", + "\n", + "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." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the pumps:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.pump.add(\n", + " Node(13, Point(4.0, -1.0)),\n", + " [pump.Static(flow_rate=[0.5 / 3600])],\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "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 minimum crest is specified at or just below the setpoint. In this example, we've chosen for the Outlet where we specify the minimum crest level 5 cm below the setpoint. For example: the Outlet of polder basin 1 (node 4) is specified with a minimum crest level of 1.95 meter." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the outlets:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up outlet\n", + "model.outlet.add(\n", + " Node(2, Point(0.0, -1.0)),\n", + " [outlet.Static(flow_rate=[2 * 0.5 / 3600], min_crest_level=[0.0])],\n", + ")\n", + "model.outlet.add(\n", + " Node(5, Point(0.0, -3.0)),\n", + " [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[1.95])],\n", + ")\n", + "model.outlet.add(\n", + " Node(7, Point(1.0, -4.0)),\n", + " [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[1.45])],\n", + ")\n", + "model.outlet.add(\n", + " Node(9, Point(3.0, -4.0)),\n", + " [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[0.95])],\n", + ")\n", + "model.outlet.add(\n", + " Node(11, Point(4.0, -3.0)),\n", + " [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[0.45])],\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "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." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the PID control:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pid_control_data = {\n", + " \"listen_node_type\": \"Basin\",\n", + " \"proportional\": [0.05],\n", + " \"integral\": [0.00],\n", + " \"derivative\": [0.0],\n", + "}\n", + "model.pid_control.add(\n", + " Node(3, Point(-1.0, -1.0)),\n", + " [pid_control.Static(listen_node_id=[4], target=[2.0], **pid_control_data)],\n", + ")\n", + "model.pid_control.add(\n", + " Node(14, Point(-1.0, -3.0)),\n", + " [pid_control.Static(listen_node_id=[6], target=[1.5], **pid_control_data)],\n", + ")\n", + "model.pid_control.add(\n", + " Node(15, Point(1.0, -3.0)),\n", + " [pid_control.Static(listen_node_id=[8], target=[1.0], **pid_control_data)],\n", + ")\n", + "model.pid_control.add(\n", + " Node(16, Point(3.0, -3.0)),\n", + " [pid_control.Static(listen_node_id=[10], target=[0.5], **pid_control_data)],\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the edges:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.edge.add(model.basin[1], model.outlet[2])\n", + "model.edge.add(model.pid_control[3], model.outlet[2])\n", + "model.edge.add(model.outlet[2], model.basin[4])\n", + "model.edge.add(model.basin[4], model.outlet[5])\n", + "model.edge.add(model.outlet[5], model.basin[6])\n", + "model.edge.add(model.basin[6], model.outlet[7])\n", + "model.edge.add(model.outlet[7], model.basin[8])\n", + "model.edge.add(model.basin[8], model.outlet[9])\n", + "model.edge.add(model.outlet[9], model.basin[10])\n", + "model.edge.add(model.basin[10], model.outlet[11])\n", + "model.edge.add(model.outlet[11], model.basin[12])\n", + "model.edge.add(model.basin[12], model.pump[13])\n", + "model.edge.add(model.pump[13], model.basin[1])\n", + "model.edge.add(model.pid_control[14], model.outlet[5])\n", + "model.edge.add(model.pid_control[15], model.outlet[7])\n", + "model.edge.add(model.pid_control[16], model.outlet[9])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To plot the model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Write the model to a TOML file and run it in the Julia." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "datadir = Path(\"data\")\n", + "model.write(datadir / \"local_pidcontrolled_cascade/ribasim.toml\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# | include: false\n", + "from subprocess import run\n", + "\n", + "run(\n", + " [\n", + " \"julia\",\n", + " \"--project=../../core\",\n", + " \"--eval\",\n", + " f'using Ribasim; Ribasim.main(\"{datadir.as_posix()}/local_pidcontrolled_cascade/ribasim.toml\")',\n", + " ],\n", + " check=True,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After running the model, read back the result to plot the flow of each polder basin." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "datadir_flow = datadir / \"local_pidcontrolled_cascade/results/flow.arrow\"\n", + "df_flow = pd.read_feather(datadir_flow)\n", + "df_flow[\"edge\"] = list(zip(df_flow.from_node_id, df_flow.to_node_id))\n", + "df_flow[\"flow_m3d\"] = df_flow.flow_rate * 86400\n", + "\n", + "df_pivot = df_flow.pivot_table(index=\"time\", columns=\"edge\", values=\"flow_m3d\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df_input = df_pivot.loc[:, [(1, 2), (13, 1)]]\n", + "df_input.plot(ylim=[-1.0, 20.0])\n", + "df_weirs = df_pivot.loc[:, [(4, 5), (6, 7), (8, 9), (10, 11)]]\n", + "df_weirs.plot(ylim=[-1.0, 15.0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " Below graph shows the vertical flux on each basin." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "datadir_basin = datadir / \"local_pidcontrolled_cascade/results/basin.arrow\"\n", + "df_basin = pd.read_feather(datadir_basin)\n", + "df_basin[\"vertical_flux\"] = (\n", + " df_basin[\"precipitation\"]\n", + " - df_basin[\"evaporation\"]\n", + " + df_basin[\"drainage\"]\n", + " + df_basin[\"infiltration\"]\n", + ")\n", + "df_basin_wide = df_basin.pivot_table(\n", + " index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\", \"vertical_flux\"]\n", + ")\n", + "df_basin_wide[\"vertical_flux\"].plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df_basin_wide[\"level\"].plot()" + ] } ], "metadata": { diff --git a/docs/python/image/cascade-polder.jpg b/docs/python/image/cascade-polder.jpg new file mode 100644 index 000000000..7fff88d49 Binary files /dev/null and b/docs/python/image/cascade-polder.jpg differ diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index e4dd9dbad..ceeb1fad2 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -34,6 +34,7 @@ pump_discrete_control_model, tabulated_rating_curve_control_model, ) +from ribasim_testmodels.doc_example import local_pidcontrolled_cascade_model from ribasim_testmodels.dutch_waterways import dutch_waterways_model from ribasim_testmodels.equations import ( linear_resistance_model, @@ -79,6 +80,7 @@ "level_setpoint_with_minmax_model", "linear_resistance_demand_model", "linear_resistance_model", + "local_pidcontrolled_cascade_model", "looped_subnetwork_model", "main_network_with_subnetworks_model", "manning_resistance_model", diff --git a/python/ribasim_testmodels/ribasim_testmodels/doc_example.py b/python/ribasim_testmodels/ribasim_testmodels/doc_example.py new file mode 100644 index 000000000..ea6a562e1 --- /dev/null +++ b/python/ribasim_testmodels/ribasim_testmodels/doc_example.py @@ -0,0 +1,178 @@ +import numpy as np +import pandas as pd +from ribasim import Model, Node +from ribasim.nodes import ( + basin, + outlet, + pid_control, + pump, +) +from shapely.geometry import Point + + +def local_pidcontrolled_cascade_model(): + """Demonstrating model for the cascade polder project from our partner""" + + model = Model(starttime="2020-01-01", endtime="2021-01-01", crs="EPSG:28992") + + # Set up 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 + drainage3 = drainage.copy() + drainage3[210:240] = 17 * 2.87e-8 + drainage4 = drainage.copy() + drainage4[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, + urban_runoff=0.0, + ), + ] + + basin_time4 = [ + basin.Time( + time=pd.date_range(model.starttime, model.endtime), + drainage=drainage4, + potential_evaporation=evaporation, + infiltration=0.0, + precipitation=precipitation, + urban_runoff=0.0, + ), + ] + basin_time3 = [ + basin.Time( + time=pd.date_range(model.starttime, model.endtime), + drainage=drainage3, + potential_evaporation=evaporation, + infiltration=0.0, + precipitation=precipitation, + urban_runoff=0.0, + ), + ] + + 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, + urban_runoff=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], + ) + + # Set up pid control + pid_control_data = { + "listen_node_type": "Basin", + "proportional": [0.1], + "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)], + ) + + # Set up pump + model.pump.add( + Node(13, Point(4.0, -1.0)), + [pump.Static(flow_rate=[0.5 / 3600])], + ) + + # Set up outlet + model.outlet.add( + Node(2, Point(0.0, -1.0)), + [outlet.Static(flow_rate=[4 * 0.5 / 3600], min_crest_level=[0.0])], + ) + model.outlet.add( + Node(5, Point(0.0, -3.0)), + [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[1.95])], + ) + model.outlet.add( + Node(7, Point(1.0, -4.0)), + [outlet.Static(flow_rate=[4 * 0.5 / 3600], min_crest_level=[1.45])], + ) + model.outlet.add( + Node(9, Point(3.0, -4.0)), + [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[0.95])], + ) + model.outlet.add( + Node(11, Point(4.0, -3.0)), + [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[0.45])], + ) + + 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]) + + return model