diff --git a/concept/allocation.html b/concept/allocation.html index 11d35cd8a..fbb5167c9 100644 --- a/concept/allocation.html +++ b/concept/allocation.html @@ -588,7 +588,7 @@

4.4 Example

The following is an example of an optimization problem for the example shown here:

-
+
Code 
using Ribasim
@@ -611,39 +611,39 @@ 
println(p.allocation.allocation_models[1].problem)
-
Min F[(Basin #5, UserDemand #6)]² + F[(Basin #2, UserDemand #3)]² + F[(Basin #12, UserDemand #13)]²
+
Min F[(Basin #5, UserDemand #6)]² + F[(Basin #12, UserDemand #13)]² + F[(Basin #2, UserDemand #3)]²
 Subject to
  source[(FlowBoundary #1, Basin #2)] : F[(FlowBoundary #1, Basin #2)] ≤ 172800
- source_user[UserDemand #6] : F[(UserDemand #6, Basin #5)] ≤ 0
- source_user[UserDemand #13] : F[(UserDemand #13, Terminal #10)] ≤ 0
  source_user[UserDemand #3] : F[(UserDemand #3, Basin #2)] ≤ 0
- fractional_flow[(TabulatedRatingCurve #7, FractionalFlow #8)] : -0.6 F[(Basin #5, TabulatedRatingCurve #7)] + F[(TabulatedRatingCurve #7, FractionalFlow #8)] ≤ 0
+ source_user[UserDemand #13] : F[(UserDemand #13, Terminal #10)] ≤ 0
+ source_user[UserDemand #6] : F[(UserDemand #6, Basin #5)] ≤ 0
+ fractional_flow[(TabulatedRatingCurve #7, FractionalFlow #8)] : F[(TabulatedRatingCurve #7, FractionalFlow #8)] - 0.6 F[(Basin #5, TabulatedRatingCurve #7)] ≤ 0
  fractional_flow[(TabulatedRatingCurve #7, FractionalFlow #9)] : -0.4 F[(Basin #5, TabulatedRatingCurve #7)] + F[(TabulatedRatingCurve #7, FractionalFlow #9)] ≤ 0
  flow_buffer_outflow[TabulatedRatingCurve #7] : F_flow_buffer_out[TabulatedRatingCurve #7] ≤ 0
- flow_conservation[Basin #2] : F[(FlowBoundary #1, Basin #2)] - F[(Basin #2, UserDemand #3)] - F[(Basin #2, LinearResistance #4)] + F[(LinearResistance #4, Basin #2)] + F[(UserDemand #3, Basin #2)] = 0
- flow_conservation[FractionalFlow #9] : -F[(FractionalFlow #9, Basin #12)] + F[(TabulatedRatingCurve #7, FractionalFlow #9)] = 0
- flow_conservation[Basin #12] : F[(FractionalFlow #9, Basin #12)] - F[(Basin #12, UserDemand #13)] = 0
- flow_conservation[Terminal #10] : F[(FractionalFlow #8, Terminal #10)] + F[(UserDemand #13, Terminal #10)] = 0
- flow_conservation[Basin #5] : -F[(Basin #5, TabulatedRatingCurve #7)] - F[(Basin #5, UserDemand #6)] + F[(UserDemand #6, Basin #5)] + F[(LinearResistance #4, Basin #5)] - F[(Basin #5, LinearResistance #4)] = 0
+ flow_conservation[Basin #2] : -F[(Basin #2, LinearResistance #4)] + F[(LinearResistance #4, Basin #2)] + F[(FlowBoundary #1, Basin #2)] + F[(UserDemand #3, Basin #2)] - F[(Basin #2, UserDemand #3)] = 0
+ flow_conservation[Basin #5] : -F[(Basin #5, UserDemand #6)] + F[(UserDemand #6, Basin #5)] + F[(LinearResistance #4, Basin #5)] - F[(Basin #5, LinearResistance #4)] - F[(Basin #5, TabulatedRatingCurve #7)] = 0
+ flow_conservation[Terminal #10] : F[(UserDemand #13, Terminal #10)] + F[(FractionalFlow #8, Terminal #10)] = 0
  flow_conservation[FractionalFlow #8] : -F[(FractionalFlow #8, Terminal #10)] + F[(TabulatedRatingCurve #7, FractionalFlow #8)] = 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[TabulatedRatingCurve #7] : F[(Basin #5, TabulatedRatingCurve #7)] - F[(TabulatedRatingCurve #7, FractionalFlow #9)] - F[(TabulatedRatingCurve #7, FractionalFlow #8)] - F_flow_buffer_in[TabulatedRatingCurve #7] + F_flow_buffer_out[TabulatedRatingCurve #7] = 0
- F[(Basin #5, TabulatedRatingCurve #7)] ≥ 0
- F[(FlowBoundary #1, Basin #2)] ≥ 0
+ flow_conservation[FractionalFlow #9] : -F[(FractionalFlow #9, Basin #12)] + F[(TabulatedRatingCurve #7, FractionalFlow #9)] = 0
+ flow_conservation[TabulatedRatingCurve #7] : -F[(TabulatedRatingCurve #7, FractionalFlow #8)] + F[(Basin #5, TabulatedRatingCurve #7)] - F[(TabulatedRatingCurve #7, FractionalFlow #9)] - F_flow_buffer_in[TabulatedRatingCurve #7] + F_flow_buffer_out[TabulatedRatingCurve #7] = 0
+ flow_conservation[Basin #12] : -F[(Basin #12, UserDemand #13)] + F[(FractionalFlow #9, Basin #12)] = 0
  F[(Basin #5, UserDemand #6)] ≥ 0
- F[(FractionalFlow #8, Terminal #10)] ≥ 0
- F[(Basin #2, UserDemand #3)] ≥ 0
  F[(Basin #2, LinearResistance #4)] ≥ 0
  F[(LinearResistance #4, Basin #2)] ≥ 0
- F[(UserDemand #3, Basin #2)] ≥ 0
- F[(UserDemand #6, Basin #5)] ≥ 0
- F[(FractionalFlow #9, Basin #12)] ≥ 0
- F[(TabulatedRatingCurve #7, FractionalFlow #9)] ≥ 0
+ F[(FlowBoundary #1, Basin #2)] ≥ 0
+ F[(UserDemand #13, Terminal #10)] ≥ 0
  F[(Basin #12, UserDemand #13)] ≥ 0
+ F[(UserDemand #6, Basin #5)] ≥ 0
  F[(LinearResistance #4, Basin #5)] ≥ 0
  F[(Basin #5, LinearResistance #4)] ≥ 0
+ F[(FractionalFlow #8, Terminal #10)] ≥ 0
  F[(TabulatedRatingCurve #7, FractionalFlow #8)] ≥ 0
- F[(UserDemand #13, Terminal #10)] ≥ 0
+ F[(FractionalFlow #9, Basin #12)] ≥ 0
+ F[(UserDemand #3, Basin #2)] ≥ 0
+ F[(Basin #5, TabulatedRatingCurve #7)] ≥ 0
+ F[(TabulatedRatingCurve #7, FractionalFlow #9)] ≥ 0
+ F[(Basin #2, UserDemand #3)] ≥ 0
  F_flow_buffer_in[TabulatedRatingCurve #7] ≥ 0
  F_flow_buffer_out[TabulatedRatingCurve #7] ≥ 0
 

diff --git a/concept/modelconcept.html b/concept/modelconcept.html
index 4f7b4a1fa..feddcc0dc 100644
--- a/concept/modelconcept.html
+++ b/concept/modelconcept.html
@@ -319,9 +319,9 @@ 

 
\[
 \frac{\mathrm{d}S}{\mathrm{d}t} = Q_{in} - Q_{out}
 \]

-
We can split out the fluxes into separate terms, such as precipitation \(P\), evapotranspiration \(ET\) and runoff \(R\). For now other fluxes are combined into \(Q_{rest}\). If we define all fluxes entering our reservoir as positive, and those leaving the system as negative, all fluxes can be summed up.

+
We can split out the fluxes into separate terms, such as precipitation \(P\) and evapotranspiration \(ET\). For now other fluxes are combined into \(Q_{rest}\). If we define all fluxes entering our reservoir as positive, and those leaving the system as negative, all fluxes can be summed up.

 
\[
-\frac{\mathrm{d}S}{\mathrm{d}t} = R + P + ET + Q_{rest}
+\frac{\mathrm{d}S}{\mathrm{d}t} = P + ET + Q_{rest}
 \]

 
 

diff --git a/guide/examples.html b/guide/examples.html
index 668dfde4f..73688521d 100644
--- a/guide/examples.html
+++ b/guide/examples.html
@@ -327,15 +327,14 @@ 
1 Basic model wit
         potential_evaporation=evaporation,
         infiltration=0.0,
         precipitation=precipitation,
-        urban_runoff=0.0,
-    ),
-    basin.State(level=[1.4]),
-]
-
-model.basin.add(Node(1, Point(0.0, 0.0)), basin_data)
-model.basin.add(Node(3, Point(2.0, 0.0)), basin_data)
-model.basin.add(Node(6, Point(3.0, 2.0)), basin_data)
-model.basin.add(Node(9, Point(5.0, 0.0)), basin_data)
+ ), + basin.State(level=[1.4]), +] + +model.basin.add(Node(1, Point(0.0, 0.0)), basin_data) +model.basin.add(Node(3, Point(2.0, 0.0)), basin_data) +model.basin.add(Node(6, Point(3.0, 2.0)), basin_data) +model.basin.add(Node(9, Point(5.0, 0.0)), basin_data)

Setup linear resistance:

@@ -1004,23 +1003,21 @@

5 Model with allo potential_evaporation=[0.0, 0.0], infiltration=[0.0, 0.0], precipitation=[1e-6, 0.0], - urban_runoff=[0.0, 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], - urban_runoff=[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], + ), + ], +)

Setup the flow boundary:

@@ -1150,66 +1147,62 @@

6 Guidance of mod 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=drainage_4, - 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=drainage_3, - 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, - ), - ], + ), +] + +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(4, Point(0.0, -2.0)), - [basin.State(level=[1.5]), polder_profile, *basin_time], + Node(6, Point(0.0, -4.0)), + [basin.State(level=[1.0]), polder_profile, *basin_time], ) model.basin.add( - Node(6, Point(0.0, -4.0)), - [basin.State(level=[1.0]), polder_profile, *basin_time], + Node(8, Point(2.0, -4.0)), + [basin.State(level=[1.5]), polder_profile, *basin_time3], ) model.basin.add( - Node(8, Point(2.0, -4.0)), - [basin.State(level=[1.5]), polder_profile, *basin_time3], + Node(10, Point(4.0, -4.0)), + [basin.State(level=[1.3]), polder_profile, *basin_time4], ) 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], -)

+ Node(12, Point(4.0, -2.0)), + [basin.State(level=[0.1]), polder_profile, *basin_time], +)

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.

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.

diff --git a/guide/examples_files/figure-html/cell-19-output-1.png b/guide/examples_files/figure-html/cell-19-output-1.png index 00e092f92..f0011f379 100644 Binary files a/guide/examples_files/figure-html/cell-19-output-1.png and b/guide/examples_files/figure-html/cell-19-output-1.png differ diff --git a/guide/examples_files/figure-html/cell-20-output-1.png b/guide/examples_files/figure-html/cell-20-output-1.png index e6c78c23e..d1f1d0e0d 100644 Binary files a/guide/examples_files/figure-html/cell-20-output-1.png and b/guide/examples_files/figure-html/cell-20-output-1.png differ diff --git a/guide/examples_files/figure-html/cell-59-output-1.png b/guide/examples_files/figure-html/cell-59-output-1.png index 54561750f..f4ea4178f 100644 Binary files a/guide/examples_files/figure-html/cell-59-output-1.png and b/guide/examples_files/figure-html/cell-59-output-1.png differ diff --git a/guide/examples_files/figure-html/cell-60-output-2.png b/guide/examples_files/figure-html/cell-60-output-2.png index 239bc10a4..80efb3eb9 100644 Binary files a/guide/examples_files/figure-html/cell-60-output-2.png and b/guide/examples_files/figure-html/cell-60-output-2.png differ diff --git a/guide/examples_files/figure-html/cell-71-output-1.png b/guide/examples_files/figure-html/cell-71-output-1.png index 01e7f9117..f13573534 100644 Binary files a/guide/examples_files/figure-html/cell-71-output-1.png and b/guide/examples_files/figure-html/cell-71-output-1.png differ diff --git a/reference/node/basin.html b/reference/node/basin.html index 7e06620b5..f8958a468 100644 --- a/reference/node/basin.html +++ b/reference/node/basin.html @@ -441,12 +441,6 @@

\(m^3 s^{-1}\) non-negative - -urban_runoff -Float64 -\(m^3 s^{-1}\) -non-negative -

Note that if variables are not set in the static table, default values are used when possible. These are generally zero, e.g. no precipitation, no inflow. If it is not possible to have a reasonable and safe default, a value must be provided in the static table.

@@ -832,7 +826,7 @@

Here \(p > 0\) is the threshold value which determines the interval \([0,p]\) of the smooth transition between \(0\) and \(1\), see the plot below.

-
+
Code 
import numpy as np
diff --git a/reference/node/discrete-control.html b/reference/node/discrete-control.html
index a6fa553df..b7ffa9baa 100644
--- a/reference/node/discrete-control.html
+++ b/reference/node/discrete-control.html
@@ -383,7 +383,7 @@ 
DiscreteControl

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

-

+

 

 Code
 
using Ribasim
diff --git a/reference/test-models.html b/reference/test-models.html
index 7fd8d61a5..66c3ed27c 100644
--- a/reference/test-models.html
+++ b/reference/test-models.html
@@ -343,7 +343,7 @@ 
Test models

 
 
 
Ribasim developers use the following models in their testbench and in order to test new features.

-

+

 

 Code
 
import ribasim_testmodels
diff --git a/reference/validation.html b/reference/validation.html
index c7cf6d9c4..d845c7c4a 100644
--- a/reference/validation.html
+++ b/reference/validation.html
@@ -353,7 +353,7 @@ 
Validation

 

 
1 Connectivity

 
In the table below, each column shows which node types are allowed to be downstream (or ‘down-control’) of the node type at the top of the column.

-

+

 

 Code
 
using Ribasim
@@ -703,7 +703,7 @@ 
1 Connectivity
 
2 Neighbor amounts

 
The table below shows for each node type between which bounds the amount of in- and outneighbors must be, for both flow and control edges.

-

+

 

 Code
 
flow_in_min = Vector{String}()
diff --git a/search.json b/search.json
index 29aebc5d8..01c12cd46 100644
--- a/search.json
+++ b/search.json
@@ -48,7 +48,7 @@
     "href": "guide/examples.html",
     "title": "Examples",
     "section": "",
-    "text": "1 Basic model with static forcing\n\nimport shutil\nfrom pathlib import Path\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nfrom ribasim import Allocation, Model, Node\nfrom ribasim.nodes import (\n    basin,\n    discrete_control,\n    flow_boundary,\n    fractional_flow,\n    level_boundary,\n    level_demand,\n    linear_resistance,\n    manning_resistance,\n    outlet,\n    pid_control,\n    pump,\n    tabulated_rating_curve,\n    user_demand,\n)\nfrom shapely.geometry import Point\n\n\ndatadir = Path(\"data\")\nshutil.rmtree(datadir, ignore_errors=True)\n\n\nmodel = Model(starttime=\"2020-01-01\", endtime=\"2021-01-01\", crs=\"EPSG:4326\")\n\nSetup the basins:\n\ntime = pd.date_range(model.starttime, model.endtime)\nday_of_year = time.day_of_year.to_numpy()\nseconds_per_day = 24 * 60 * 60\nevaporation = (\n    (-1.0 * np.cos(day_of_year / 365.0 * 2 * np.pi) + 1.0) * 0.0025 / seconds_per_day\n)\nrng = np.random.default_rng(seed=0)\nprecipitation = (\n    rng.lognormal(mean=-1.0, sigma=1.7, size=time.size) * 0.001 / seconds_per_day\n)\n\n# Convert steady forcing to m/s\n# 2 mm/d precipitation, 1 mm/d evaporation\n\nbasin_data = [\n    basin.Profile(area=[0.01, 1000.0], level=[0.0, 1.0]),\n    basin.Time(\n        time=pd.date_range(model.starttime, model.endtime),\n        drainage=0.0,\n        potential_evaporation=evaporation,\n        infiltration=0.0,\n        precipitation=precipitation,\n        urban_runoff=0.0,\n    ),\n    basin.State(level=[1.4]),\n]\n\nmodel.basin.add(Node(1, Point(0.0, 0.0)), basin_data)\nmodel.basin.add(Node(3, Point(2.0, 0.0)), basin_data)\nmodel.basin.add(Node(6, Point(3.0, 2.0)), basin_data)\nmodel.basin.add(Node(9, Point(5.0, 0.0)), basin_data)\n\nSetup linear resistance:\n\nmodel.linear_resistance.add(\n    Node(10, Point(6.0, 0.0)),\n    [linear_resistance.Static(resistance=[5e3])],\n)\nmodel.linear_resistance.add(\n    Node(12, Point(2.0, 1.0)),\n    [linear_resistance.Static(resistance=[3600.0 * 24.0 / 100.0])],\n)\n\nSetup Manning resistance:\n\nmodel.manning_resistance.add(\n    Node(2, Point(1.0, 0.0)),\n    [\n        manning_resistance.Static(\n            length=[900], manning_n=[0.04], profile_width=[6.0], profile_slope=[3.0]\n        )\n    ],\n)\n\nSet up a rating curve node:\n\nmodel.tabulated_rating_curve.add(\n    Node(4, Point(3.0, 0.0)),\n    [tabulated_rating_curve.Static(level=[0.0, 1.0], flow_rate=[0.0, 10 / 86400])],\n)\n\nSetup fractional flows:\n\nmodel.fractional_flow.add(\n    Node(5, Point(3.0, 1.0)), [fractional_flow.Static(fraction=[0.3])]\n)\nmodel.fractional_flow.add(\n    Node(8, Point(4.0, 0.0)), [fractional_flow.Static(fraction=[0.6])]\n)\nmodel.fractional_flow.add(\n    Node(13, Point(3.0, -1.0)),\n    [fractional_flow.Static(fraction=[0.1])],\n)\n\nSetup pump:\n\nmodel.pump.add(Node(7, Point(4.0, 1.0)), [pump.Static(flow_rate=[0.5 / 3600])])\n\nSetup level boundary:\n\nmodel.level_boundary.add(\n    Node(11, Point(2.0, 2.0)), [level_boundary.Static(level=[0.5])]\n)\nmodel.level_boundary.add(\n    Node(17, Point(6.0, 1.0)), [level_boundary.Static(level=[1.5])]\n)\n\nSetup flow boundary:\n\nmodel.flow_boundary.add(\n    Node(15, Point(3.0, 3.0)), [flow_boundary.Static(flow_rate=[1e-4])]\n)\nmodel.flow_boundary.add(\n    Node(16, Point(0.0, 1.0)), [flow_boundary.Static(flow_rate=[1e-4])]\n)\n\nSetup terminal:\n\nmodel.terminal.add(Node(14, Point(3.0, -2.0)))\n\nSetup the edges:\n\nmodel.edge.add(model.basin[1], model.manning_resistance[2])\nmodel.edge.add(model.manning_resistance[2], model.basin[3])\nmodel.edge.add(model.basin[3], model.tabulated_rating_curve[4])\nmodel.edge.add(model.tabulated_rating_curve[4], model.fractional_flow[5])\nmodel.edge.add(model.tabulated_rating_curve[4], model.fractional_flow[8])\nmodel.edge.add(model.fractional_flow[5], model.basin[6])\nmodel.edge.add(model.basin[6], model.pump[7])\nmodel.edge.add(model.fractional_flow[8], model.basin[9])\nmodel.edge.add(model.pump[7], model.basin[9])\nmodel.edge.add(model.basin[9], model.linear_resistance[10])\nmodel.edge.add(model.level_boundary[11], model.linear_resistance[12])\nmodel.edge.add(model.linear_resistance[12], model.basin[3])\nmodel.edge.add(model.tabulated_rating_curve[4], model.fractional_flow[13])\nmodel.edge.add(model.fractional_flow[13], model.terminal[14])\nmodel.edge.add(model.flow_boundary[15], model.basin[6])\nmodel.edge.add(model.flow_boundary[16], model.basin[1])\nmodel.edge.add(model.linear_resistance[10], model.level_boundary[17])\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n\n\n\n\n\n\n\nWrite the model to a TOML and GeoPackage:\n\ntoml_path = datadir / \"basic/ribasim.toml\"\nmodel.write(toml_path)\n\nPosixPath('data/basic/ribasim.toml')\n\n\nNow run the model. From Python you can run it with:\nimport subprocess\nsubprocess.run([cli_path, toml_path], check=True)\nOr similarly you can from the terminal with:\nribasim basic/ribasim.toml\nAfter running the model, read back the results:\n\ndf_basin = pd.read_feather(datadir / \"basic/results/basin.arrow\")\ndf_basin_wide = df_basin.pivot_table(\n    index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\ndf_basin_wide[\"level\"].plot()\n\n\n\n\n\n\n\n\n\ndf_flow = pd.read_feather(datadir / \"basic/results/flow.arrow\")\ndf_flow[\"edge\"] = list(zip(df_flow.from_node_id, df_flow.to_node_id))\ndf_flow[\"flow_m3d\"] = df_flow.flow_rate * 86400\nax = df_flow.pivot_table(index=\"time\", columns=\"edge\", values=\"flow_m3d\").plot()\nax.legend(bbox_to_anchor=(1.3, 1), title=\"Edge\")\n\n\n\n\n\n\n\n\n\n\n2 Model with discrete control\nThe 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.\nSetup the basins:\n\nmodel = Model(starttime=\"2020-01-01\", endtime=\"2021-01-01\", crs=\"EPSG:4326\")\n\n\nmodel.basin.add(\n    Node(1, Point(0.0, 0.0)),\n    [\n        basin.Profile(area=[1000.0, 1000.0], level=[0.0, 1.0]),\n        basin.State(level=[20.0]),\n        basin.Time(time=[\"2020-01-01\", \"2020-07-01\"], precipitation=[0.0, 3e-6]),\n    ],\n)\n\nSetup the discrete control:\n\nmodel.discrete_control.add(\n    Node(7, Point(1.0, 0.0)),\n    [\n        discrete_control.Variable(\n            compound_variable_id=1,\n            listen_node_id=1,\n            listen_node_type=[\"Basin\"],\n            variable=[\"level\"],\n        ),\n        discrete_control.Condition(\n            compound_variable_id=1,\n            # min, max\n            greater_than=[5.0, 15.0],\n        ),\n        discrete_control.Logic(\n            truth_state=[\"FF\", \"TF\", \"TT\"],\n            control_state=[\"in\", \"none\", \"out\"],\n        ),\n    ],\n)\n\nThe above control logic can be summarized as follows:\n\nIf the level is above the maximum, activate the control state “out”;\nIf the level is below the minimum, active the control state “in”;\nOtherwise activate the control state “none”.\n\nSetup the pump:\n\nmodel.pump.add(\n    Node(2, Point(1.0, 1.0)),\n    [pump.Static(control_state=[\"none\", \"in\", \"out\"], flow_rate=[0.0, 2e-3, 0.0])],\n)\nmodel.pump.add(\n    Node(3, Point(1.0, -1.0)),\n    [pump.Static(control_state=[\"none\", \"in\", \"out\"], flow_rate=[0.0, 0.0, 2e-3])],\n)\n\nThe pump data defines the following:\n\n\n\nControl state\nPump #2 flow rate (m/s)\nPump #3 flow rate (m/s)\n\n\n\n\n“none”\n0.0\n0.0\n\n\n“in”\n2e-3\n0.0\n\n\n“out”\n0.0\n2e-3\n\n\n\nSetup the level boundary:\n\nmodel.level_boundary.add(\n    Node(4, Point(2.0, 0.0)), [level_boundary.Static(level=[10.0])]\n)\n\nSetup the rating curve:\n\nmodel.tabulated_rating_curve.add(\n    Node(5, Point(-1.0, 0.0)),\n    [tabulated_rating_curve.Static(level=[2.0, 15.0], flow_rate=[0.0, 2e-3])],\n)\n\nSetup the terminal:\n\nmodel.terminal.add(Node(6, Point(-2.0, 0.0)))\n\nSetup edges:\n\nmodel.edge.add(model.basin[1], model.pump[3])\nmodel.edge.add(model.pump[3], model.level_boundary[4])\nmodel.edge.add(model.level_boundary[4], model.pump[2])\nmodel.edge.add(model.pump[2], model.basin[1])\nmodel.edge.add(model.basin[1], model.tabulated_rating_curve[5])\nmodel.edge.add(model.tabulated_rating_curve[5], model.terminal[6])\nmodel.edge.add(model.discrete_control[7], model.pump[2])\nmodel.edge.add(model.discrete_control[7], model.pump[3])\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n\n\n\n\n\n\n\nListen edges are plotted with a dashed line since they are not present in the “Edge / static” schema but only in the “Control / condition” schema.\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"level_range/ribasim.toml\")\n\nPosixPath('data/level_range/ribasim.toml')\n\n\nNow run the model with ribasim level_range/ribasim.toml. After running the model, read back the results:\n\ndf_basin = pd.read_feather(datadir / \"level_range/results/basin.arrow\")\ndf_basin_wide = df_basin.pivot_table(\n    index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\n\nax = df_basin_wide[\"level\"].plot()\n\ngreater_than = model.discrete_control.condition.df.greater_than\n\nax.hlines(\n    greater_than,\n    df_basin.time[0],\n    df_basin.time.max(),\n    lw=1,\n    ls=\"--\",\n    color=\"k\",\n)\n\nax.set_yticks(greater_than, [\"min\", \"max\"])\nax.set_ylabel(\"level\")\nplt.show()\n\n\n\n\n\n\n\n\nWe see that in January the level of the basin is too high and thus water is pumped out until the maximum level of the desired range is reached. Then until May water flows out of the basin freely through the tabulated rating curve until the minimum level is reached. From this point until the start of July water is pumped into the basin in short bursts to stay within the desired range. At the start of July rain starts falling on the basin, which causes the basin level to rise until the maximum level. From this point onward water is pumped out of the basin in short bursts to stay within the desired range.\n\n\n3 Model with PID control\nSet up the model:\n\nmodel = Model(starttime=\"2020-01-01\", endtime=\"2020-12-01\", crs=\"EPSG:4326\")\n\nSetup the basins:\n\nmodel.basin.add(\n    Node(2, Point(1.0, 0.0)),\n    [basin.Profile(area=[1000.0, 1000.0], level=[0.0, 1.0]), basin.State(level=[6.0])],\n)\n\nSetup the pump:\n\nmodel.pump.add(\n    Node(3, Point(2.0, 0.5)),\n    [pump.Static(flow_rate=[0.0])],  # Will be overwritten by PID controller\n)\n\nSetup the outlet:\n\nmodel.outlet.add(\n    Node(6, Point(2.0, -0.5)),\n    [outlet.Static(flow_rate=[0.0])],  # Will be overwritten by PID controller\n)\n\nSetup flow boundary:\n\nmodel.flow_boundary.add(\n    Node(1, Point(0.0, 0.0)),\n    [flow_boundary.Static(flow_rate=[1e-3])],\n)\n\nSetup level boundary:\n\nmodel.level_boundary.add(\n    Node(4, Point(3.0, 0.0)),\n    [level_boundary.Static(level=[5.0])],\n)\n\nSetup PID control:\n\nfor node, proportional, integral in [\n    (Node(5, Point(1.5, 1.0)), -1e-3, -1e-7),\n    (Node(7, Point(1.5, -1.0)), 1e-3, 1e-7),\n]:\n    pid_control_data = [\n        pid_control.Time(\n            time=[\n                \"2020-01-01\",\n                \"2020-05-01\",\n                \"2020-07-01\",\n                \"2020-12-01\",\n            ],\n            listen_node_id=2,\n            listen_node_type=\"Basin\",\n            target=[5.0, 5.0, 7.5, 7.5],\n            proportional=proportional,\n            integral=integral,\n            derivative=0.0,\n        )\n    ]\n    model.pid_control.add(node, pid_control_data)\n\nNote 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.\nSetup the edges:\n\nmodel.edge.add(model.flow_boundary[1], model.basin[2])\nmodel.edge.add(model.basin[2], model.pump[3])\nmodel.edge.add(model.pump[3], model.level_boundary[4])\nmodel.edge.add(model.level_boundary[4], model.outlet[6])\nmodel.edge.add(model.outlet[6], model.basin[2])\nmodel.edge.add(model.pid_control[5], model.pump[3])\nmodel.edge.add(model.pid_control[7], model.outlet[6])\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n\n\n\n\n\n\n\nWrite the model to a TOML and GeoPackage:\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"pid_control/ribasim.toml\")\n\nPosixPath('data/pid_control/ribasim.toml')\n\n\nNow run the model with ribasim pid_control/ribasim.toml. After running the model, read back the results:\n\nfrom matplotlib.dates import date2num\n\ndf_basin = pd.read_feather(datadir / \"pid_control/results/basin.arrow\")\ndf_basin_wide = df_basin.pivot_table(\n    index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\nax = df_basin_wide[\"level\"].plot()\nax.set_ylabel(\"level [m]\")\n\n# Plot target level\nlevel_demands = model.pid_control.time.df.target.to_numpy()[:4]\ntimes = date2num(model.pid_control.time.df.time)[:4]\nax.plot(times, level_demands, color=\"k\", ls=\":\", label=\"target level\")\npass\n\n\n\n\n\n\n\n\n\n\n4 Model with allocation (user demand)\nSetup a model:\n\nmodel = Model(starttime=\"2020-01-01\", endtime=\"2020-01-20\", crs=\"EPSG:4326\")\n\nSetup the basins:\n\nbasin_data = [\n    basin.Profile(area=[300_000.0, 300_000.0], level=[0.0, 1.0]),\n    basin.State(level=[1.0]),\n]\n\nmodel.basin.add(\n    Node(2, Point(1.0, 0.0), subnetwork_id=1),\n    basin_data,\n)\nmodel.basin.add(\n    Node(5, Point(3.0, 0.0), subnetwork_id=1),\n    basin_data,\n)\nmodel.basin.add(\n    Node(12, Point(4.5, 1.0), subnetwork_id=1),\n    basin_data,\n)\n\nSetup the flow boundary:\n\nmodel.flow_boundary.add(\n    Node(1, Point(0.0, 0.0), subnetwork_id=1), [flow_boundary.Static(flow_rate=[2.0])]\n)\n\nSetup the linear resistance:\n\nmodel.linear_resistance.add(\n    Node(4, Point(2.0, 0.0), subnetwork_id=1),\n    [linear_resistance.Static(resistance=[0.06])],\n)\n\nSetup the tabulated rating curve:\n\nmodel.tabulated_rating_curve.add(\n    Node(7, Point(4.0, 0.0), subnetwork_id=1),\n    [tabulated_rating_curve.Static(level=[0.0, 0.5, 1.0], flow_rate=[0.0, 0.0, 2.0])],\n)\n\nSetup the fractional flow:\n\nmodel.fractional_flow.add(\n    Node(8, Point(4.5, 0.0), subnetwork_id=1),\n    [fractional_flow.Static(fraction=[0.6, 0.9], control_state=[\"divert\", \"close\"])],\n)\nmodel.fractional_flow.add(\n    Node(9, Point(4.5, 0.5), subnetwork_id=1),\n    [fractional_flow.Static(fraction=[0.4, 0.1], control_state=[\"divert\", \"close\"])],\n)\n\nSetup the terminal:\n\nmodel.terminal.add(Node(10, Point(5.0, 0.0), subnetwork_id=1))\n\nSetup the discrete control:\n\nmodel.discrete_control.add(\n    Node(11, Point(4.5, 0.25), subnetwork_id=1),\n    [\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            truth_state=[\"T\", \"F\"], control_state=[\"divert\", \"close\"]\n        ),\n    ],\n)\n\nSetup the users:\n\nmodel.user_demand.add(\n    Node(6, Point(3.0, 1.0), subnetwork_id=1),\n    [\n        user_demand.Static(\n            demand=[1.5], return_factor=[0.0], min_level=[-1.0], priority=[1]\n        )\n    ],\n)\nmodel.user_demand.add(\n    Node(13, Point(5.0, 1.0), subnetwork_id=1),\n    [\n        user_demand.Static(\n            demand=[1.0], return_factor=[0.0], min_level=[-1.0], priority=[3]\n        )\n    ],\n)\nmodel.user_demand.add(\n    Node(3, Point(1.0, 1.0), subnetwork_id=1),\n    [\n        user_demand.Time(\n            demand=[0.0, 1.0, 1.2, 1.2],\n            return_factor=[0.0, 0.0, 0.0, 0.0],\n            min_level=[-1.0, -1.0, -1.0, -1.0],\n            priority=[1, 1, 2, 2],\n            time=2 * [\"2020-01-01\", \"2020-01-20\"],\n        )\n    ],\n)\n\nSetup the allocation:\n\nmodel.allocation = Allocation(use_allocation=True, timestep=86400)\n\nSetup the edges:\n\nmodel.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=1)\nmodel.edge.add(model.basin[2], model.user_demand[3])\nmodel.edge.add(model.basin[2], model.linear_resistance[4])\nmodel.edge.add(model.linear_resistance[4], model.basin[5])\nmodel.edge.add(model.basin[5], model.user_demand[6])\nmodel.edge.add(model.basin[5], model.tabulated_rating_curve[7])\nmodel.edge.add(model.tabulated_rating_curve[7], model.fractional_flow[8])\nmodel.edge.add(model.user_demand[3], model.basin[2])\nmodel.edge.add(model.user_demand[6], model.basin[5])\nmodel.edge.add(model.tabulated_rating_curve[7], model.fractional_flow[9])\nmodel.edge.add(model.fractional_flow[8], model.terminal[10])\nmodel.edge.add(model.fractional_flow[9], model.basin[12])\nmodel.edge.add(model.basin[12], model.user_demand[13])\nmodel.edge.add(model.user_demand[13], model.terminal[10])\nmodel.edge.add(model.discrete_control[11], model.fractional_flow[8])\nmodel.edge.add(model.discrete_control[11], model.fractional_flow[9])\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n\n\n\n\n\n\n\nWrite the model to a TOML and GeoPackage:\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"allocation_example/ribasim.toml\")\n\nPosixPath('data/allocation_example/ribasim.toml')\n\n\nNow run the model with ribasim allocation_example/ribasim.toml. After running the model, read back the results:\n\nimport matplotlib.ticker as plticker\n\ndf_allocation = pd.read_feather(datadir / \"allocation_example/results/allocation.arrow\")\ndf_allocation_wide = df_allocation.pivot_table(\n    index=\"time\",\n    columns=[\"node_type\", \"node_id\", \"priority\"],\n    values=[\"demand\", \"allocated\", \"realized\"],\n)\ndf_allocation_wide = df_allocation_wide.loc[:, (df_allocation_wide != 0).any(axis=0)]\n\nfig, axs = plt.subplots(1, 3, figsize=(8, 5))\n\ndf_allocation_wide[\"demand\"].plot(ax=axs[0], ls=\":\")\ndf_allocation_wide[\"allocated\"].plot(ax=axs[1], ls=\"--\")\ndf_allocation_wide.xs(1, level=\"priority\", axis=1)[\"realized\"].plot(\n    ax=axs[2], color=[\"C0\", \"C2\", \"C3\"]\n)\n\nfig.tight_layout()\nloc = plticker.MultipleLocator(2)\n\naxs[0].set_ylabel(\"level [m]\")\n\nfor ax, title in zip(axs, [\"Demand\", \"Allocated\", \"Abstracted\"]):\n    ax.set_title(title)\n    ax.set_ylim(0.0, 1.6)\n    ax.xaxis.set_major_locator(loc)\n\n\n\n\n\n\n\n\nSome things to note about this plot:\n\nThe realized flow at the start time is not correct, as there is no realized flow yet. The given value is how much of its total demand a user can abstract in the physical layer.\nNo flow was allocated to UserDemand 13 so that is not plotted.\nAbstraction is accumulated over all priorities per user.\n\n\ndf_basin = pd.read_feather(datadir / \"allocation_example/results/basin.arrow\")\ndf_basin_wide = df_basin.pivot_table(\n    index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\n\nax = df_basin_wide[\"level\"].plot()\nax.set_title(\"Basin levels\")\nax.set_ylabel(\"level [m]\")\n\nText(0, 0.5, 'level [m]')\n\n\n\n\n\n\n\n\n\n\n\n5 Model with allocation (basin supply/demand)\nSetup a model:\n\nmodel = Model(starttime=\"2020-01-01\", endtime=\"2020-02-01\", crs=\"EPSG:4326\")\n\nSetup the basins:\n\nbasin_data = [\n    basin.Profile(area=[1e3, 1e3], level=[0.0, 1.0]),\n    basin.State(level=[0.5]),\n]\nmodel.basin.add(\n    Node(2, Point(1.0, 0.0), subnetwork_id=2),\n    [\n        *basin_data,\n        basin.Time(\n            time=[\"2020-01-01\", \"2020-01-16\"],\n            drainage=[0.0, 0.0],\n            potential_evaporation=[0.0, 0.0],\n            infiltration=[0.0, 0.0],\n            precipitation=[1e-6, 0.0],\n            urban_runoff=[0.0, 0.0],\n        ),\n    ],\n)\nmodel.basin.add(\n    Node(5, Point(2.0, -1.0), subnetwork_id=2),\n    [\n        *basin_data,\n        basin.Static(\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\nSetup the flow boundary:\n\nmodel.flow_boundary.add(\n    Node(1, Point(0.0, 0.0), subnetwork_id=2), [flow_boundary.Static(flow_rate=[1e-3])]\n)\n\nSetup level demand:\n\nmodel.level_demand.add(\n    Node(4, Point(1.0, -1.0), subnetwork_id=2),\n    [level_demand.Static(priority=[1], min_level=[1.0], max_level=[1.5])],\n)\n\nSetup the users:\n\nmodel.user_demand.add(\n    Node(3, Point(2.0, 0.0), subnetwork_id=2),\n    [\n        user_demand.Static(\n            priority=[2], demand=[1.5e-3], return_factor=[0.2], min_level=[0.2]\n        )\n    ],\n)\n\nSetup the allocation:\n\nmodel.allocation = Allocation(use_allocation=True, timestep=1e5)\n\nSetup the edges:\n\nmodel.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2)\nmodel.edge.add(model.basin[2], model.user_demand[3])\nmodel.edge.add(model.level_demand[4], model.basin[2])\nmodel.edge.add(model.user_demand[3], model.basin[5])\nmodel.edge.add(model.level_demand[4], model.basin[5])\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n\n\n\n\n\n\n\nWrite the model to a TOML and GeoPackage:\n\nmodel.write(datadir / \"level_demand/ribasim.toml\")\n\nPosixPath('data/level_demand/ribasim.toml')\n\n\nNow run the model with ribasim level_demand/ribasim.toml. After running the model, read back the results:\n\ndf_basin = pd.read_feather(datadir / \"level_demand/results/basin.arrow\")\ndf_basin = df_basin[df_basin.node_id == 2]\ndf_basin_wide = df_basin.pivot_table(\n    index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\nax = df_basin_wide[\"level\"].plot(ylabel=\"level [m]\")\n\n\n\n\n\n\n\n\nIn the plot above, the line denotes the level of Basin #2 over time. The Basin level is a piecewise linear function of time, with several stages explained below.\nConstants:\n\n\\(d\\): UserDemand #3 demand,\n\\(\\phi\\): Basin #2 precipitation rate,\n\\(q\\): LevelBoundary flow.\n\nStages:\n\nIn the first stage the Basin takes precedence so the UserDemand doesn’t abstract, hence the net change of Basin #2 is \\(q + \\phi\\);\nIn the second stage (and following stages) the Basin no longer has a positive demand, since precipitation provides enough water to get the Basin to its target level. The FlowBoundary flow gets fully allocated to the UserDemand, hence the net change of Basin #2 is \\(\\phi\\);\nIn the third stage the Basin enters its surplus stage, even though initially the level is below the maximum level. This is because the simulation anticipates that the current precipitation is going to bring the Basin level over its maximum level. The net change of Basin #2 is now \\(q + \\phi - d\\);\nAt the start of the fourth stage the precipitation stops, and so the UserDemand partly uses surplus water from the Basin to fulfill its demand. The net change of Basin #2 becomes \\(q - d\\).\nIn the final stage the Basin is in a dynamical equilibrium, since the Basin has no supply so the user abstracts precisely the flow from the LevelBoundary.\n\n\n\n6 Guidance of modelling a cascade of polder basins\nSituation 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.\n\nThe 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\n\nmodel = Model(starttime=\"2020-01-01\", endtime=\"2021-01-01\", crs=\"EPSG:28992\")\n\nAll the polder basins are exposed to time varying forcings (precipitation, evaporation, drainage, infiltration) to mimic situations of water excess and water shortage.\nIn 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.\nModel 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.\nSetting up the basins:\n\ntime = pd.date_range(model.starttime, model.endtime)\nday_of_year = time.day_of_year.to_numpy()\n\nprecipitation = np.zeros(day_of_year.size)\nprecipitation[0:90] = 1.72e-8\nprecipitation[330:366] = 1.72e-8\n\nevaporation = np.zeros(day_of_year.size)\nevaporation[130:270] = 2.87e-8\n\ndrainage = np.zeros(day_of_year.size)\ndrainage[120:270] = 0.4 * 2.87e-8\ndrainage_3 = drainage.copy()\ndrainage_3[210:240] = 17 * 2.87e-8\ndrainage_4 = drainage.copy()\ndrainage_4[160:240] = 13 * 2.87e-8\n\ninfiltration = np.zeros(day_of_year.size)\ninfiltration[0:90] = 5e-8\n\npolder_profile = basin.Profile(area=[100, 100], level=[0.0, 3.0])\n\nbasin_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\nbasin_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]\nbasin_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\nmodel.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)\nmodel.basin.add(\n    Node(4, Point(0.0, -2.0)),\n    [basin.State(level=[1.5]), polder_profile, *basin_time],\n)\nmodel.basin.add(\n    Node(6, Point(0.0, -4.0)),\n    [basin.State(level=[1.0]), polder_profile, *basin_time],\n)\nmodel.basin.add(\n    Node(8, Point(2.0, -4.0)),\n    [basin.State(level=[1.5]), polder_profile, *basin_time3],\n)\nmodel.basin.add(\n    Node(10, Point(4.0, -4.0)),\n    [basin.State(level=[1.3]), polder_profile, *basin_time4],\n)\nmodel.basin.add(\n    Node(12, Point(4.0, -2.0)),\n    [basin.State(level=[0.1]), polder_profile, *basin_time],\n)\n\nAfter 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.\nSince 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.\nSetup the pumps:\n\nmodel.pump.add(\n    Node(13, Point(4.0, -1.0)),\n    [pump.Static(flow_rate=[0.5 / 3600])],\n)\n\nAccording 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.\nSetup the outlets:\n\n# Set up outlet\nmodel.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)\nmodel.outlet.add(\n    Node(5, Point(0.0, -3.0)),\n    [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[1.95])],\n)\nmodel.outlet.add(\n    Node(7, Point(1.0, -4.0)),\n    [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[1.45])],\n)\nmodel.outlet.add(\n    Node(9, Point(3.0, -4.0)),\n    [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[0.95])],\n)\nmodel.outlet.add(\n    Node(11, Point(4.0, -3.0)),\n    [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[0.45])],\n)\n\nWhen 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.\nSetup the PID control:\n\npid_control_data = {\n    \"listen_node_type\": \"Basin\",\n    \"proportional\": [0.05],\n    \"integral\": [0.00],\n    \"derivative\": [0.0],\n}\nmodel.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)\nmodel.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)\nmodel.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)\nmodel.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)\n\nSetup the edges:\n\nmodel.edge.add(model.basin[1], model.outlet[2])\nmodel.edge.add(model.pid_control[3], model.outlet[2])\nmodel.edge.add(model.outlet[2], model.basin[4])\nmodel.edge.add(model.basin[4], model.outlet[5])\nmodel.edge.add(model.outlet[5], model.basin[6])\nmodel.edge.add(model.basin[6], model.outlet[7])\nmodel.edge.add(model.outlet[7], model.basin[8])\nmodel.edge.add(model.basin[8], model.outlet[9])\nmodel.edge.add(model.outlet[9], model.basin[10])\nmodel.edge.add(model.basin[10], model.outlet[11])\nmodel.edge.add(model.outlet[11], model.basin[12])\nmodel.edge.add(model.basin[12], model.pump[13])\nmodel.edge.add(model.pump[13], model.basin[1])\nmodel.edge.add(model.pid_control[14], model.outlet[5])\nmodel.edge.add(model.pid_control[15], model.outlet[7])\nmodel.edge.add(model.pid_control[16], model.outlet[9])\n\nTo plot the model\n\nmodel.plot()\n\n\n\n\n\n\n\n\nWrite the model to a TOML file and run it in the Julia.\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"local_pidcontrolled_cascade/ribasim.toml\")\n\nPosixPath('data/local_pidcontrolled_cascade/ribasim.toml')\n\n\nAfter running the model, read back the result to plot the flow of each polder basin.\n\ndatadir_flow = datadir / \"local_pidcontrolled_cascade/results/flow.arrow\"\ndf_flow = pd.read_feather(datadir_flow)\ndf_flow[\"edge\"] = list(zip(df_flow.from_node_id, df_flow.to_node_id))\ndf_flow[\"flow_m3d\"] = df_flow.flow_rate * 86400\n\ndf_pivot = df_flow.pivot_table(index=\"time\", columns=\"edge\", values=\"flow_m3d\")\n\nBelow graphs show the flow exchanged with the mainsystem (i.e. the inlet and the pump), and the flow of weirs inbetween the polder basins.\n\ndf_input = df_pivot.loc[:, [(1, 2), (13, 1)]]\ndf_input.plot(ylim=[-1.0, 20.0])\ndf_weirs = df_pivot.loc[:, [(4, 5), (6, 7), (8, 9), (10, 11)]]\ndf_weirs.plot(ylim=[-1.0, 15.0])\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nBelow graph shows the vertical flux on each basin.\n\ndatadir_basin = datadir / \"local_pidcontrolled_cascade/results/basin.arrow\"\ndf_basin = pd.read_feather(datadir_basin)\ndf_basin[\"vertical_flux\"] = (\n    df_basin[\"precipitation\"]\n    - df_basin[\"evaporation\"]\n    + df_basin[\"drainage\"]\n    + df_basin[\"infiltration\"]\n)\ndf_basin_wide = df_basin.pivot_table(\n    index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\", \"vertical_flux\"]\n)\ndf_basin_wide[\"vertical_flux\"].plot()\n\n\n\n\n\n\n\n\nIn 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.\n\ndf_basin_wide[\"level\"].plot()",
+    "text": "1 Basic model with static forcing\n\nimport shutil\nfrom pathlib import Path\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nfrom ribasim import Allocation, Model, Node\nfrom ribasim.nodes import (\n    basin,\n    discrete_control,\n    flow_boundary,\n    fractional_flow,\n    level_boundary,\n    level_demand,\n    linear_resistance,\n    manning_resistance,\n    outlet,\n    pid_control,\n    pump,\n    tabulated_rating_curve,\n    user_demand,\n)\nfrom shapely.geometry import Point\n\n\ndatadir = Path(\"data\")\nshutil.rmtree(datadir, ignore_errors=True)\n\n\nmodel = Model(starttime=\"2020-01-01\", endtime=\"2021-01-01\", crs=\"EPSG:4326\")\n\nSetup the basins:\n\ntime = pd.date_range(model.starttime, model.endtime)\nday_of_year = time.day_of_year.to_numpy()\nseconds_per_day = 24 * 60 * 60\nevaporation = (\n    (-1.0 * np.cos(day_of_year / 365.0 * 2 * np.pi) + 1.0) * 0.0025 / seconds_per_day\n)\nrng = np.random.default_rng(seed=0)\nprecipitation = (\n    rng.lognormal(mean=-1.0, sigma=1.7, size=time.size) * 0.001 / seconds_per_day\n)\n\n# Convert steady forcing to m/s\n# 2 mm/d precipitation, 1 mm/d evaporation\n\nbasin_data = [\n    basin.Profile(area=[0.01, 1000.0], level=[0.0, 1.0]),\n    basin.Time(\n        time=pd.date_range(model.starttime, model.endtime),\n        drainage=0.0,\n        potential_evaporation=evaporation,\n        infiltration=0.0,\n        precipitation=precipitation,\n    ),\n    basin.State(level=[1.4]),\n]\n\nmodel.basin.add(Node(1, Point(0.0, 0.0)), basin_data)\nmodel.basin.add(Node(3, Point(2.0, 0.0)), basin_data)\nmodel.basin.add(Node(6, Point(3.0, 2.0)), basin_data)\nmodel.basin.add(Node(9, Point(5.0, 0.0)), basin_data)\n\nSetup linear resistance:\n\nmodel.linear_resistance.add(\n    Node(10, Point(6.0, 0.0)),\n    [linear_resistance.Static(resistance=[5e3])],\n)\nmodel.linear_resistance.add(\n    Node(12, Point(2.0, 1.0)),\n    [linear_resistance.Static(resistance=[3600.0 * 24.0 / 100.0])],\n)\n\nSetup Manning resistance:\n\nmodel.manning_resistance.add(\n    Node(2, Point(1.0, 0.0)),\n    [\n        manning_resistance.Static(\n            length=[900], manning_n=[0.04], profile_width=[6.0], profile_slope=[3.0]\n        )\n    ],\n)\n\nSet up a rating curve node:\n\nmodel.tabulated_rating_curve.add(\n    Node(4, Point(3.0, 0.0)),\n    [tabulated_rating_curve.Static(level=[0.0, 1.0], flow_rate=[0.0, 10 / 86400])],\n)\n\nSetup fractional flows:\n\nmodel.fractional_flow.add(\n    Node(5, Point(3.0, 1.0)), [fractional_flow.Static(fraction=[0.3])]\n)\nmodel.fractional_flow.add(\n    Node(8, Point(4.0, 0.0)), [fractional_flow.Static(fraction=[0.6])]\n)\nmodel.fractional_flow.add(\n    Node(13, Point(3.0, -1.0)),\n    [fractional_flow.Static(fraction=[0.1])],\n)\n\nSetup pump:\n\nmodel.pump.add(Node(7, Point(4.0, 1.0)), [pump.Static(flow_rate=[0.5 / 3600])])\n\nSetup level boundary:\n\nmodel.level_boundary.add(\n    Node(11, Point(2.0, 2.0)), [level_boundary.Static(level=[0.5])]\n)\nmodel.level_boundary.add(\n    Node(17, Point(6.0, 1.0)), [level_boundary.Static(level=[1.5])]\n)\n\nSetup flow boundary:\n\nmodel.flow_boundary.add(\n    Node(15, Point(3.0, 3.0)), [flow_boundary.Static(flow_rate=[1e-4])]\n)\nmodel.flow_boundary.add(\n    Node(16, Point(0.0, 1.0)), [flow_boundary.Static(flow_rate=[1e-4])]\n)\n\nSetup terminal:\n\nmodel.terminal.add(Node(14, Point(3.0, -2.0)))\n\nSetup the edges:\n\nmodel.edge.add(model.basin[1], model.manning_resistance[2])\nmodel.edge.add(model.manning_resistance[2], model.basin[3])\nmodel.edge.add(model.basin[3], model.tabulated_rating_curve[4])\nmodel.edge.add(model.tabulated_rating_curve[4], model.fractional_flow[5])\nmodel.edge.add(model.tabulated_rating_curve[4], model.fractional_flow[8])\nmodel.edge.add(model.fractional_flow[5], model.basin[6])\nmodel.edge.add(model.basin[6], model.pump[7])\nmodel.edge.add(model.fractional_flow[8], model.basin[9])\nmodel.edge.add(model.pump[7], model.basin[9])\nmodel.edge.add(model.basin[9], model.linear_resistance[10])\nmodel.edge.add(model.level_boundary[11], model.linear_resistance[12])\nmodel.edge.add(model.linear_resistance[12], model.basin[3])\nmodel.edge.add(model.tabulated_rating_curve[4], model.fractional_flow[13])\nmodel.edge.add(model.fractional_flow[13], model.terminal[14])\nmodel.edge.add(model.flow_boundary[15], model.basin[6])\nmodel.edge.add(model.flow_boundary[16], model.basin[1])\nmodel.edge.add(model.linear_resistance[10], model.level_boundary[17])\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n\n\n\n\n\n\n\nWrite the model to a TOML and GeoPackage:\n\ntoml_path = datadir / \"basic/ribasim.toml\"\nmodel.write(toml_path)\n\nPosixPath('data/basic/ribasim.toml')\n\n\nNow run the model. From Python you can run it with:\nimport subprocess\nsubprocess.run([cli_path, toml_path], check=True)\nOr similarly you can from the terminal with:\nribasim basic/ribasim.toml\nAfter running the model, read back the results:\n\ndf_basin = pd.read_feather(datadir / \"basic/results/basin.arrow\")\ndf_basin_wide = df_basin.pivot_table(\n    index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\ndf_basin_wide[\"level\"].plot()\n\n\n\n\n\n\n\n\n\ndf_flow = pd.read_feather(datadir / \"basic/results/flow.arrow\")\ndf_flow[\"edge\"] = list(zip(df_flow.from_node_id, df_flow.to_node_id))\ndf_flow[\"flow_m3d\"] = df_flow.flow_rate * 86400\nax = df_flow.pivot_table(index=\"time\", columns=\"edge\", values=\"flow_m3d\").plot()\nax.legend(bbox_to_anchor=(1.3, 1), title=\"Edge\")\n\n\n\n\n\n\n\n\n\n\n2 Model with discrete control\nThe 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.\nSetup the basins:\n\nmodel = Model(starttime=\"2020-01-01\", endtime=\"2021-01-01\", crs=\"EPSG:4326\")\n\n\nmodel.basin.add(\n    Node(1, Point(0.0, 0.0)),\n    [\n        basin.Profile(area=[1000.0, 1000.0], level=[0.0, 1.0]),\n        basin.State(level=[20.0]),\n        basin.Time(time=[\"2020-01-01\", \"2020-07-01\"], precipitation=[0.0, 3e-6]),\n    ],\n)\n\nSetup the discrete control:\n\nmodel.discrete_control.add(\n    Node(7, Point(1.0, 0.0)),\n    [\n        discrete_control.Variable(\n            compound_variable_id=1,\n            listen_node_id=1,\n            listen_node_type=[\"Basin\"],\n            variable=[\"level\"],\n        ),\n        discrete_control.Condition(\n            compound_variable_id=1,\n            # min, max\n            greater_than=[5.0, 15.0],\n        ),\n        discrete_control.Logic(\n            truth_state=[\"FF\", \"TF\", \"TT\"],\n            control_state=[\"in\", \"none\", \"out\"],\n        ),\n    ],\n)\n\nThe above control logic can be summarized as follows:\n\nIf the level is above the maximum, activate the control state “out”;\nIf the level is below the minimum, active the control state “in”;\nOtherwise activate the control state “none”.\n\nSetup the pump:\n\nmodel.pump.add(\n    Node(2, Point(1.0, 1.0)),\n    [pump.Static(control_state=[\"none\", \"in\", \"out\"], flow_rate=[0.0, 2e-3, 0.0])],\n)\nmodel.pump.add(\n    Node(3, Point(1.0, -1.0)),\n    [pump.Static(control_state=[\"none\", \"in\", \"out\"], flow_rate=[0.0, 0.0, 2e-3])],\n)\n\nThe pump data defines the following:\n\n\n\nControl state\nPump #2 flow rate (m/s)\nPump #3 flow rate (m/s)\n\n\n\n\n“none”\n0.0\n0.0\n\n\n“in”\n2e-3\n0.0\n\n\n“out”\n0.0\n2e-3\n\n\n\nSetup the level boundary:\n\nmodel.level_boundary.add(\n    Node(4, Point(2.0, 0.0)), [level_boundary.Static(level=[10.0])]\n)\n\nSetup the rating curve:\n\nmodel.tabulated_rating_curve.add(\n    Node(5, Point(-1.0, 0.0)),\n    [tabulated_rating_curve.Static(level=[2.0, 15.0], flow_rate=[0.0, 2e-3])],\n)\n\nSetup the terminal:\n\nmodel.terminal.add(Node(6, Point(-2.0, 0.0)))\n\nSetup edges:\n\nmodel.edge.add(model.basin[1], model.pump[3])\nmodel.edge.add(model.pump[3], model.level_boundary[4])\nmodel.edge.add(model.level_boundary[4], model.pump[2])\nmodel.edge.add(model.pump[2], model.basin[1])\nmodel.edge.add(model.basin[1], model.tabulated_rating_curve[5])\nmodel.edge.add(model.tabulated_rating_curve[5], model.terminal[6])\nmodel.edge.add(model.discrete_control[7], model.pump[2])\nmodel.edge.add(model.discrete_control[7], model.pump[3])\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n\n\n\n\n\n\n\nListen edges are plotted with a dashed line since they are not present in the “Edge / static” schema but only in the “Control / condition” schema.\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"level_range/ribasim.toml\")\n\nPosixPath('data/level_range/ribasim.toml')\n\n\nNow run the model with ribasim level_range/ribasim.toml. After running the model, read back the results:\n\ndf_basin = pd.read_feather(datadir / \"level_range/results/basin.arrow\")\ndf_basin_wide = df_basin.pivot_table(\n    index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\n\nax = df_basin_wide[\"level\"].plot()\n\ngreater_than = model.discrete_control.condition.df.greater_than\n\nax.hlines(\n    greater_than,\n    df_basin.time[0],\n    df_basin.time.max(),\n    lw=1,\n    ls=\"--\",\n    color=\"k\",\n)\n\nax.set_yticks(greater_than, [\"min\", \"max\"])\nax.set_ylabel(\"level\")\nplt.show()\n\n\n\n\n\n\n\n\nWe see that in January the level of the basin is too high and thus water is pumped out until the maximum level of the desired range is reached. Then until May water flows out of the basin freely through the tabulated rating curve until the minimum level is reached. From this point until the start of July water is pumped into the basin in short bursts to stay within the desired range. At the start of July rain starts falling on the basin, which causes the basin level to rise until the maximum level. From this point onward water is pumped out of the basin in short bursts to stay within the desired range.\n\n\n3 Model with PID control\nSet up the model:\n\nmodel = Model(starttime=\"2020-01-01\", endtime=\"2020-12-01\", crs=\"EPSG:4326\")\n\nSetup the basins:\n\nmodel.basin.add(\n    Node(2, Point(1.0, 0.0)),\n    [basin.Profile(area=[1000.0, 1000.0], level=[0.0, 1.0]), basin.State(level=[6.0])],\n)\n\nSetup the pump:\n\nmodel.pump.add(\n    Node(3, Point(2.0, 0.5)),\n    [pump.Static(flow_rate=[0.0])],  # Will be overwritten by PID controller\n)\n\nSetup the outlet:\n\nmodel.outlet.add(\n    Node(6, Point(2.0, -0.5)),\n    [outlet.Static(flow_rate=[0.0])],  # Will be overwritten by PID controller\n)\n\nSetup flow boundary:\n\nmodel.flow_boundary.add(\n    Node(1, Point(0.0, 0.0)),\n    [flow_boundary.Static(flow_rate=[1e-3])],\n)\n\nSetup level boundary:\n\nmodel.level_boundary.add(\n    Node(4, Point(3.0, 0.0)),\n    [level_boundary.Static(level=[5.0])],\n)\n\nSetup PID control:\n\nfor node, proportional, integral in [\n    (Node(5, Point(1.5, 1.0)), -1e-3, -1e-7),\n    (Node(7, Point(1.5, -1.0)), 1e-3, 1e-7),\n]:\n    pid_control_data = [\n        pid_control.Time(\n            time=[\n                \"2020-01-01\",\n                \"2020-05-01\",\n                \"2020-07-01\",\n                \"2020-12-01\",\n            ],\n            listen_node_id=2,\n            listen_node_type=\"Basin\",\n            target=[5.0, 5.0, 7.5, 7.5],\n            proportional=proportional,\n            integral=integral,\n            derivative=0.0,\n        )\n    ]\n    model.pid_control.add(node, pid_control_data)\n\nNote 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.\nSetup the edges:\n\nmodel.edge.add(model.flow_boundary[1], model.basin[2])\nmodel.edge.add(model.basin[2], model.pump[3])\nmodel.edge.add(model.pump[3], model.level_boundary[4])\nmodel.edge.add(model.level_boundary[4], model.outlet[6])\nmodel.edge.add(model.outlet[6], model.basin[2])\nmodel.edge.add(model.pid_control[5], model.pump[3])\nmodel.edge.add(model.pid_control[7], model.outlet[6])\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n\n\n\n\n\n\n\nWrite the model to a TOML and GeoPackage:\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"pid_control/ribasim.toml\")\n\nPosixPath('data/pid_control/ribasim.toml')\n\n\nNow run the model with ribasim pid_control/ribasim.toml. After running the model, read back the results:\n\nfrom matplotlib.dates import date2num\n\ndf_basin = pd.read_feather(datadir / \"pid_control/results/basin.arrow\")\ndf_basin_wide = df_basin.pivot_table(\n    index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\nax = df_basin_wide[\"level\"].plot()\nax.set_ylabel(\"level [m]\")\n\n# Plot target level\nlevel_demands = model.pid_control.time.df.target.to_numpy()[:4]\ntimes = date2num(model.pid_control.time.df.time)[:4]\nax.plot(times, level_demands, color=\"k\", ls=\":\", label=\"target level\")\npass\n\n\n\n\n\n\n\n\n\n\n4 Model with allocation (user demand)\nSetup a model:\n\nmodel = Model(starttime=\"2020-01-01\", endtime=\"2020-01-20\", crs=\"EPSG:4326\")\n\nSetup the basins:\n\nbasin_data = [\n    basin.Profile(area=[300_000.0, 300_000.0], level=[0.0, 1.0]),\n    basin.State(level=[1.0]),\n]\n\nmodel.basin.add(\n    Node(2, Point(1.0, 0.0), subnetwork_id=1),\n    basin_data,\n)\nmodel.basin.add(\n    Node(5, Point(3.0, 0.0), subnetwork_id=1),\n    basin_data,\n)\nmodel.basin.add(\n    Node(12, Point(4.5, 1.0), subnetwork_id=1),\n    basin_data,\n)\n\nSetup the flow boundary:\n\nmodel.flow_boundary.add(\n    Node(1, Point(0.0, 0.0), subnetwork_id=1), [flow_boundary.Static(flow_rate=[2.0])]\n)\n\nSetup the linear resistance:\n\nmodel.linear_resistance.add(\n    Node(4, Point(2.0, 0.0), subnetwork_id=1),\n    [linear_resistance.Static(resistance=[0.06])],\n)\n\nSetup the tabulated rating curve:\n\nmodel.tabulated_rating_curve.add(\n    Node(7, Point(4.0, 0.0), subnetwork_id=1),\n    [tabulated_rating_curve.Static(level=[0.0, 0.5, 1.0], flow_rate=[0.0, 0.0, 2.0])],\n)\n\nSetup the fractional flow:\n\nmodel.fractional_flow.add(\n    Node(8, Point(4.5, 0.0), subnetwork_id=1),\n    [fractional_flow.Static(fraction=[0.6, 0.9], control_state=[\"divert\", \"close\"])],\n)\nmodel.fractional_flow.add(\n    Node(9, Point(4.5, 0.5), subnetwork_id=1),\n    [fractional_flow.Static(fraction=[0.4, 0.1], control_state=[\"divert\", \"close\"])],\n)\n\nSetup the terminal:\n\nmodel.terminal.add(Node(10, Point(5.0, 0.0), subnetwork_id=1))\n\nSetup the discrete control:\n\nmodel.discrete_control.add(\n    Node(11, Point(4.5, 0.25), subnetwork_id=1),\n    [\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            truth_state=[\"T\", \"F\"], control_state=[\"divert\", \"close\"]\n        ),\n    ],\n)\n\nSetup the users:\n\nmodel.user_demand.add(\n    Node(6, Point(3.0, 1.0), subnetwork_id=1),\n    [\n        user_demand.Static(\n            demand=[1.5], return_factor=[0.0], min_level=[-1.0], priority=[1]\n        )\n    ],\n)\nmodel.user_demand.add(\n    Node(13, Point(5.0, 1.0), subnetwork_id=1),\n    [\n        user_demand.Static(\n            demand=[1.0], return_factor=[0.0], min_level=[-1.0], priority=[3]\n        )\n    ],\n)\nmodel.user_demand.add(\n    Node(3, Point(1.0, 1.0), subnetwork_id=1),\n    [\n        user_demand.Time(\n            demand=[0.0, 1.0, 1.2, 1.2],\n            return_factor=[0.0, 0.0, 0.0, 0.0],\n            min_level=[-1.0, -1.0, -1.0, -1.0],\n            priority=[1, 1, 2, 2],\n            time=2 * [\"2020-01-01\", \"2020-01-20\"],\n        )\n    ],\n)\n\nSetup the allocation:\n\nmodel.allocation = Allocation(use_allocation=True, timestep=86400)\n\nSetup the edges:\n\nmodel.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=1)\nmodel.edge.add(model.basin[2], model.user_demand[3])\nmodel.edge.add(model.basin[2], model.linear_resistance[4])\nmodel.edge.add(model.linear_resistance[4], model.basin[5])\nmodel.edge.add(model.basin[5], model.user_demand[6])\nmodel.edge.add(model.basin[5], model.tabulated_rating_curve[7])\nmodel.edge.add(model.tabulated_rating_curve[7], model.fractional_flow[8])\nmodel.edge.add(model.user_demand[3], model.basin[2])\nmodel.edge.add(model.user_demand[6], model.basin[5])\nmodel.edge.add(model.tabulated_rating_curve[7], model.fractional_flow[9])\nmodel.edge.add(model.fractional_flow[8], model.terminal[10])\nmodel.edge.add(model.fractional_flow[9], model.basin[12])\nmodel.edge.add(model.basin[12], model.user_demand[13])\nmodel.edge.add(model.user_demand[13], model.terminal[10])\nmodel.edge.add(model.discrete_control[11], model.fractional_flow[8])\nmodel.edge.add(model.discrete_control[11], model.fractional_flow[9])\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n\n\n\n\n\n\n\nWrite the model to a TOML and GeoPackage:\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"allocation_example/ribasim.toml\")\n\nPosixPath('data/allocation_example/ribasim.toml')\n\n\nNow run the model with ribasim allocation_example/ribasim.toml. After running the model, read back the results:\n\nimport matplotlib.ticker as plticker\n\ndf_allocation = pd.read_feather(datadir / \"allocation_example/results/allocation.arrow\")\ndf_allocation_wide = df_allocation.pivot_table(\n    index=\"time\",\n    columns=[\"node_type\", \"node_id\", \"priority\"],\n    values=[\"demand\", \"allocated\", \"realized\"],\n)\ndf_allocation_wide = df_allocation_wide.loc[:, (df_allocation_wide != 0).any(axis=0)]\n\nfig, axs = plt.subplots(1, 3, figsize=(8, 5))\n\ndf_allocation_wide[\"demand\"].plot(ax=axs[0], ls=\":\")\ndf_allocation_wide[\"allocated\"].plot(ax=axs[1], ls=\"--\")\ndf_allocation_wide.xs(1, level=\"priority\", axis=1)[\"realized\"].plot(\n    ax=axs[2], color=[\"C0\", \"C2\", \"C3\"]\n)\n\nfig.tight_layout()\nloc = plticker.MultipleLocator(2)\n\naxs[0].set_ylabel(\"level [m]\")\n\nfor ax, title in zip(axs, [\"Demand\", \"Allocated\", \"Abstracted\"]):\n    ax.set_title(title)\n    ax.set_ylim(0.0, 1.6)\n    ax.xaxis.set_major_locator(loc)\n\n\n\n\n\n\n\n\nSome things to note about this plot:\n\nThe realized flow at the start time is not correct, as there is no realized flow yet. The given value is how much of its total demand a user can abstract in the physical layer.\nNo flow was allocated to UserDemand 13 so that is not plotted.\nAbstraction is accumulated over all priorities per user.\n\n\ndf_basin = pd.read_feather(datadir / \"allocation_example/results/basin.arrow\")\ndf_basin_wide = df_basin.pivot_table(\n    index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\n\nax = df_basin_wide[\"level\"].plot()\nax.set_title(\"Basin levels\")\nax.set_ylabel(\"level [m]\")\n\nText(0, 0.5, 'level [m]')\n\n\n\n\n\n\n\n\n\n\n\n5 Model with allocation (basin supply/demand)\nSetup a model:\n\nmodel = Model(starttime=\"2020-01-01\", endtime=\"2020-02-01\", crs=\"EPSG:4326\")\n\nSetup the basins:\n\nbasin_data = [\n    basin.Profile(area=[1e3, 1e3], level=[0.0, 1.0]),\n    basin.State(level=[0.5]),\n]\nmodel.basin.add(\n    Node(2, Point(1.0, 0.0), subnetwork_id=2),\n    [\n        *basin_data,\n        basin.Time(\n            time=[\"2020-01-01\", \"2020-01-16\"],\n            drainage=[0.0, 0.0],\n            potential_evaporation=[0.0, 0.0],\n            infiltration=[0.0, 0.0],\n            precipitation=[1e-6, 0.0],\n        ),\n    ],\n)\nmodel.basin.add(\n    Node(5, Point(2.0, -1.0), subnetwork_id=2),\n    [\n        *basin_data,\n        basin.Static(\n            drainage=[0.0],\n            potential_evaporation=[0.0],\n            infiltration=[0.0],\n            precipitation=[0.0],\n        ),\n    ],\n)\n\nSetup the flow boundary:\n\nmodel.flow_boundary.add(\n    Node(1, Point(0.0, 0.0), subnetwork_id=2), [flow_boundary.Static(flow_rate=[1e-3])]\n)\n\nSetup level demand:\n\nmodel.level_demand.add(\n    Node(4, Point(1.0, -1.0), subnetwork_id=2),\n    [level_demand.Static(priority=[1], min_level=[1.0], max_level=[1.5])],\n)\n\nSetup the users:\n\nmodel.user_demand.add(\n    Node(3, Point(2.0, 0.0), subnetwork_id=2),\n    [\n        user_demand.Static(\n            priority=[2], demand=[1.5e-3], return_factor=[0.2], min_level=[0.2]\n        )\n    ],\n)\n\nSetup the allocation:\n\nmodel.allocation = Allocation(use_allocation=True, timestep=1e5)\n\nSetup the edges:\n\nmodel.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2)\nmodel.edge.add(model.basin[2], model.user_demand[3])\nmodel.edge.add(model.level_demand[4], model.basin[2])\nmodel.edge.add(model.user_demand[3], model.basin[5])\nmodel.edge.add(model.level_demand[4], model.basin[5])\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n\n\n\n\n\n\n\nWrite the model to a TOML and GeoPackage:\n\nmodel.write(datadir / \"level_demand/ribasim.toml\")\n\nPosixPath('data/level_demand/ribasim.toml')\n\n\nNow run the model with ribasim level_demand/ribasim.toml. After running the model, read back the results:\n\ndf_basin = pd.read_feather(datadir / \"level_demand/results/basin.arrow\")\ndf_basin = df_basin[df_basin.node_id == 2]\ndf_basin_wide = df_basin.pivot_table(\n    index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\"]\n)\nax = df_basin_wide[\"level\"].plot(ylabel=\"level [m]\")\n\n\n\n\n\n\n\n\nIn the plot above, the line denotes the level of Basin #2 over time. The Basin level is a piecewise linear function of time, with several stages explained below.\nConstants:\n\n\\(d\\): UserDemand #3 demand,\n\\(\\phi\\): Basin #2 precipitation rate,\n\\(q\\): LevelBoundary flow.\n\nStages:\n\nIn the first stage the Basin takes precedence so the UserDemand doesn’t abstract, hence the net change of Basin #2 is \\(q + \\phi\\);\nIn the second stage (and following stages) the Basin no longer has a positive demand, since precipitation provides enough water to get the Basin to its target level. The FlowBoundary flow gets fully allocated to the UserDemand, hence the net change of Basin #2 is \\(\\phi\\);\nIn the third stage the Basin enters its surplus stage, even though initially the level is below the maximum level. This is because the simulation anticipates that the current precipitation is going to bring the Basin level over its maximum level. The net change of Basin #2 is now \\(q + \\phi - d\\);\nAt the start of the fourth stage the precipitation stops, and so the UserDemand partly uses surplus water from the Basin to fulfill its demand. The net change of Basin #2 becomes \\(q - d\\).\nIn the final stage the Basin is in a dynamical equilibrium, since the Basin has no supply so the user abstracts precisely the flow from the LevelBoundary.\n\n\n\n6 Guidance of modelling a cascade of polder basins\nSituation 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.\n\nThe 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\n\nmodel = Model(starttime=\"2020-01-01\", endtime=\"2021-01-01\", crs=\"EPSG:28992\")\n\nAll the polder basins are exposed to time varying forcings (precipitation, evaporation, drainage, infiltration) to mimic situations of water excess and water shortage.\nIn 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.\nModel 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.\nSetting up the basins:\n\ntime = pd.date_range(model.starttime, model.endtime)\nday_of_year = time.day_of_year.to_numpy()\n\nprecipitation = np.zeros(day_of_year.size)\nprecipitation[0:90] = 1.72e-8\nprecipitation[330:366] = 1.72e-8\n\nevaporation = np.zeros(day_of_year.size)\nevaporation[130:270] = 2.87e-8\n\ndrainage = np.zeros(day_of_year.size)\ndrainage[120:270] = 0.4 * 2.87e-8\ndrainage_3 = drainage.copy()\ndrainage_3[210:240] = 17 * 2.87e-8\ndrainage_4 = drainage.copy()\ndrainage_4[160:240] = 13 * 2.87e-8\n\ninfiltration = np.zeros(day_of_year.size)\ninfiltration[0:90] = 5e-8\n\npolder_profile = basin.Profile(area=[100, 100], level=[0.0, 3.0])\n\nbasin_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    ),\n]\n\nbasin_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    ),\n]\nbasin_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    ),\n]\n\nmodel.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        ),\n    ],\n)\nmodel.basin.add(\n    Node(4, Point(0.0, -2.0)),\n    [basin.State(level=[1.5]), polder_profile, *basin_time],\n)\nmodel.basin.add(\n    Node(6, Point(0.0, -4.0)),\n    [basin.State(level=[1.0]), polder_profile, *basin_time],\n)\nmodel.basin.add(\n    Node(8, Point(2.0, -4.0)),\n    [basin.State(level=[1.5]), polder_profile, *basin_time3],\n)\nmodel.basin.add(\n    Node(10, Point(4.0, -4.0)),\n    [basin.State(level=[1.3]), polder_profile, *basin_time4],\n)\nmodel.basin.add(\n    Node(12, Point(4.0, -2.0)),\n    [basin.State(level=[0.1]), polder_profile, *basin_time],\n)\n\nAfter 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.\nSince 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.\nSetup the pumps:\n\nmodel.pump.add(\n    Node(13, Point(4.0, -1.0)),\n    [pump.Static(flow_rate=[0.5 / 3600])],\n)\n\nAccording 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.\nSetup the outlets:\n\n# Set up outlet\nmodel.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)\nmodel.outlet.add(\n    Node(5, Point(0.0, -3.0)),\n    [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[1.95])],\n)\nmodel.outlet.add(\n    Node(7, Point(1.0, -4.0)),\n    [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[1.45])],\n)\nmodel.outlet.add(\n    Node(9, Point(3.0, -4.0)),\n    [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[0.95])],\n)\nmodel.outlet.add(\n    Node(11, Point(4.0, -3.0)),\n    [outlet.Static(flow_rate=[0.5 / 3600], min_crest_level=[0.45])],\n)\n\nWhen 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.\nSetup the PID control:\n\npid_control_data = {\n    \"listen_node_type\": \"Basin\",\n    \"proportional\": [0.05],\n    \"integral\": [0.00],\n    \"derivative\": [0.0],\n}\nmodel.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)\nmodel.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)\nmodel.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)\nmodel.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)\n\nSetup the edges:\n\nmodel.edge.add(model.basin[1], model.outlet[2])\nmodel.edge.add(model.pid_control[3], model.outlet[2])\nmodel.edge.add(model.outlet[2], model.basin[4])\nmodel.edge.add(model.basin[4], model.outlet[5])\nmodel.edge.add(model.outlet[5], model.basin[6])\nmodel.edge.add(model.basin[6], model.outlet[7])\nmodel.edge.add(model.outlet[7], model.basin[8])\nmodel.edge.add(model.basin[8], model.outlet[9])\nmodel.edge.add(model.outlet[9], model.basin[10])\nmodel.edge.add(model.basin[10], model.outlet[11])\nmodel.edge.add(model.outlet[11], model.basin[12])\nmodel.edge.add(model.basin[12], model.pump[13])\nmodel.edge.add(model.pump[13], model.basin[1])\nmodel.edge.add(model.pid_control[14], model.outlet[5])\nmodel.edge.add(model.pid_control[15], model.outlet[7])\nmodel.edge.add(model.pid_control[16], model.outlet[9])\n\nTo plot the model\n\nmodel.plot()\n\n\n\n\n\n\n\n\nWrite the model to a TOML file and run it in the Julia.\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"local_pidcontrolled_cascade/ribasim.toml\")\n\nPosixPath('data/local_pidcontrolled_cascade/ribasim.toml')\n\n\nAfter running the model, read back the result to plot the flow of each polder basin.\n\ndatadir_flow = datadir / \"local_pidcontrolled_cascade/results/flow.arrow\"\ndf_flow = pd.read_feather(datadir_flow)\ndf_flow[\"edge\"] = list(zip(df_flow.from_node_id, df_flow.to_node_id))\ndf_flow[\"flow_m3d\"] = df_flow.flow_rate * 86400\n\ndf_pivot = df_flow.pivot_table(index=\"time\", columns=\"edge\", values=\"flow_m3d\")\n\nBelow graphs show the flow exchanged with the mainsystem (i.e. the inlet and the pump), and the flow of weirs inbetween the polder basins.\n\ndf_input = df_pivot.loc[:, [(1, 2), (13, 1)]]\ndf_input.plot(ylim=[-1.0, 20.0])\ndf_weirs = df_pivot.loc[:, [(4, 5), (6, 7), (8, 9), (10, 11)]]\ndf_weirs.plot(ylim=[-1.0, 15.0])\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\nBelow graph shows the vertical flux on each basin.\n\ndatadir_basin = datadir / \"local_pidcontrolled_cascade/results/basin.arrow\"\ndf_basin = pd.read_feather(datadir_basin)\ndf_basin[\"vertical_flux\"] = (\n    df_basin[\"precipitation\"]\n    - df_basin[\"evaporation\"]\n    + df_basin[\"drainage\"]\n    + df_basin[\"infiltration\"]\n)\ndf_basin_wide = df_basin.pivot_table(\n    index=\"time\", columns=\"node_id\", values=[\"storage\", \"level\", \"vertical_flux\"]\n)\ndf_basin_wide[\"vertical_flux\"].plot()\n\n\n\n\n\n\n\n\nIn 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.\n\ndf_basin_wide[\"level\"].plot()",
     "crumbs": [
       "How-to guides",
       "Examples"
@@ -187,7 +187,7 @@
     "href": "concept/allocation.html#example",
     "title": "Allocation",
     "section": "4.4 Example",
-    "text": "4.4 Example\nThe following is an example of an optimization problem for the example shown here:\n\n\nCode\nusing Ribasim\nusing Ribasim: NodeID\nusing SQLite\nusing ComponentArrays: ComponentVector\n\ntoml_path = normpath(@__DIR__, \"../../generated_testmodels/allocation_example/ribasim.toml\")\np = Ribasim.Model(toml_path).integrator.p\nu = ComponentVector(; storage = zeros(length(p.basin.node_id)))\n\nallocation_model = p.allocation.allocation_models[1]\nt = 0.0\npriority_idx = 1\n\nRibasim.set_flow!(p.graph, NodeID(:FlowBoundary, 1, p), NodeID(:Basin, 2, p), 1.0)\nRibasim.set_objective_priority!(allocation_model, p, u, t, priority_idx)\nRibasim.set_initial_values!(allocation_model, p, u, t)\n\nprintln(p.allocation.allocation_models[1].problem)\n\n\nMin F[(Basin #5, UserDemand #6)]² + F[(Basin #2, UserDemand #3)]² + F[(Basin #12, UserDemand #13)]²\nSubject to\n source[(FlowBoundary #1, Basin #2)] : F[(FlowBoundary #1, Basin #2)] ≤ 172800\n source_user[UserDemand #6] : F[(UserDemand #6, Basin #5)] ≤ 0\n source_user[UserDemand #13] : F[(UserDemand #13, Terminal #10)] ≤ 0\n source_user[UserDemand #3] : F[(UserDemand #3, Basin #2)] ≤ 0\n fractional_flow[(TabulatedRatingCurve #7, FractionalFlow #8)] : -0.6 F[(Basin #5, TabulatedRatingCurve #7)] + F[(TabulatedRatingCurve #7, FractionalFlow #8)] ≤ 0\n fractional_flow[(TabulatedRatingCurve #7, FractionalFlow #9)] : -0.4 F[(Basin #5, TabulatedRatingCurve #7)] + F[(TabulatedRatingCurve #7, FractionalFlow #9)] ≤ 0\n flow_buffer_outflow[TabulatedRatingCurve #7] : F_flow_buffer_out[TabulatedRatingCurve #7] ≤ 0\n flow_conservation[Basin #2] : F[(FlowBoundary #1, Basin #2)] - F[(Basin #2, UserDemand #3)] - F[(Basin #2, LinearResistance #4)] + F[(LinearResistance #4, Basin #2)] + F[(UserDemand #3, Basin #2)] = 0\n flow_conservation[FractionalFlow #9] : -F[(FractionalFlow #9, Basin #12)] + F[(TabulatedRatingCurve #7, FractionalFlow #9)] = 0\n flow_conservation[Basin #12] : F[(FractionalFlow #9, Basin #12)] - F[(Basin #12, UserDemand #13)] = 0\n flow_conservation[Terminal #10] : F[(FractionalFlow #8, Terminal #10)] + F[(UserDemand #13, Terminal #10)] = 0\n flow_conservation[Basin #5] : -F[(Basin #5, TabulatedRatingCurve #7)] - F[(Basin #5, UserDemand #6)] + F[(UserDemand #6, Basin #5)] + F[(LinearResistance #4, Basin #5)] - F[(Basin #5, LinearResistance #4)] = 0\n flow_conservation[FractionalFlow #8] : -F[(FractionalFlow #8, Terminal #10)] + F[(TabulatedRatingCurve #7, FractionalFlow #8)] = 0\n 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\n flow_conservation[TabulatedRatingCurve #7] : F[(Basin #5, TabulatedRatingCurve #7)] - F[(TabulatedRatingCurve #7, FractionalFlow #9)] - F[(TabulatedRatingCurve #7, FractionalFlow #8)] - F_flow_buffer_in[TabulatedRatingCurve #7] + F_flow_buffer_out[TabulatedRatingCurve #7] = 0\n F[(Basin #5, TabulatedRatingCurve #7)] ≥ 0\n F[(FlowBoundary #1, Basin #2)] ≥ 0\n F[(Basin #5, UserDemand #6)] ≥ 0\n F[(FractionalFlow #8, Terminal #10)] ≥ 0\n F[(Basin #2, UserDemand #3)] ≥ 0\n F[(Basin #2, LinearResistance #4)] ≥ 0\n F[(LinearResistance #4, Basin #2)] ≥ 0\n F[(UserDemand #3, Basin #2)] ≥ 0\n F[(UserDemand #6, Basin #5)] ≥ 0\n F[(FractionalFlow #9, Basin #12)] ≥ 0\n F[(TabulatedRatingCurve #7, FractionalFlow #9)] ≥ 0\n F[(Basin #12, UserDemand #13)] ≥ 0\n F[(LinearResistance #4, Basin #5)] ≥ 0\n F[(Basin #5, LinearResistance #4)] ≥ 0\n F[(TabulatedRatingCurve #7, FractionalFlow #8)] ≥ 0\n F[(UserDemand #13, Terminal #10)] ≥ 0\n F_flow_buffer_in[TabulatedRatingCurve #7] ≥ 0\n F_flow_buffer_out[TabulatedRatingCurve #7] ≥ 0",
+    "text": "4.4 Example\nThe following is an example of an optimization problem for the example shown here:\n\n\nCode\nusing Ribasim\nusing Ribasim: NodeID\nusing SQLite\nusing ComponentArrays: ComponentVector\n\ntoml_path = normpath(@__DIR__, \"../../generated_testmodels/allocation_example/ribasim.toml\")\np = Ribasim.Model(toml_path).integrator.p\nu = ComponentVector(; storage = zeros(length(p.basin.node_id)))\n\nallocation_model = p.allocation.allocation_models[1]\nt = 0.0\npriority_idx = 1\n\nRibasim.set_flow!(p.graph, NodeID(:FlowBoundary, 1, p), NodeID(:Basin, 2, p), 1.0)\nRibasim.set_objective_priority!(allocation_model, p, u, t, priority_idx)\nRibasim.set_initial_values!(allocation_model, p, u, t)\n\nprintln(p.allocation.allocation_models[1].problem)\n\n\nMin F[(Basin #5, UserDemand #6)]² + F[(Basin #12, UserDemand #13)]² + F[(Basin #2, UserDemand #3)]²\nSubject to\n source[(FlowBoundary #1, Basin #2)] : F[(FlowBoundary #1, Basin #2)] ≤ 172800\n source_user[UserDemand #3] : F[(UserDemand #3, Basin #2)] ≤ 0\n source_user[UserDemand #13] : F[(UserDemand #13, Terminal #10)] ≤ 0\n source_user[UserDemand #6] : F[(UserDemand #6, Basin #5)] ≤ 0\n fractional_flow[(TabulatedRatingCurve #7, FractionalFlow #8)] : F[(TabulatedRatingCurve #7, FractionalFlow #8)] - 0.6 F[(Basin #5, TabulatedRatingCurve #7)] ≤ 0\n fractional_flow[(TabulatedRatingCurve #7, FractionalFlow #9)] : -0.4 F[(Basin #5, TabulatedRatingCurve #7)] + F[(TabulatedRatingCurve #7, FractionalFlow #9)] ≤ 0\n flow_buffer_outflow[TabulatedRatingCurve #7] : F_flow_buffer_out[TabulatedRatingCurve #7] ≤ 0\n flow_conservation[Basin #2] : -F[(Basin #2, LinearResistance #4)] + F[(LinearResistance #4, Basin #2)] + F[(FlowBoundary #1, Basin #2)] + F[(UserDemand #3, Basin #2)] - F[(Basin #2, UserDemand #3)] = 0\n flow_conservation[Basin #5] : -F[(Basin #5, UserDemand #6)] + F[(UserDemand #6, Basin #5)] + F[(LinearResistance #4, Basin #5)] - F[(Basin #5, LinearResistance #4)] - F[(Basin #5, TabulatedRatingCurve #7)] = 0\n flow_conservation[Terminal #10] : F[(UserDemand #13, Terminal #10)] + F[(FractionalFlow #8, Terminal #10)] = 0\n flow_conservation[FractionalFlow #8] : -F[(FractionalFlow #8, Terminal #10)] + F[(TabulatedRatingCurve #7, FractionalFlow #8)] = 0\n 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\n flow_conservation[FractionalFlow #9] : -F[(FractionalFlow #9, Basin #12)] + F[(TabulatedRatingCurve #7, FractionalFlow #9)] = 0\n flow_conservation[TabulatedRatingCurve #7] : -F[(TabulatedRatingCurve #7, FractionalFlow #8)] + F[(Basin #5, TabulatedRatingCurve #7)] - F[(TabulatedRatingCurve #7, FractionalFlow #9)] - F_flow_buffer_in[TabulatedRatingCurve #7] + F_flow_buffer_out[TabulatedRatingCurve #7] = 0\n flow_conservation[Basin #12] : -F[(Basin #12, UserDemand #13)] + F[(FractionalFlow #9, Basin #12)] = 0\n F[(Basin #5, UserDemand #6)] ≥ 0\n F[(Basin #2, LinearResistance #4)] ≥ 0\n F[(LinearResistance #4, Basin #2)] ≥ 0\n F[(FlowBoundary #1, Basin #2)] ≥ 0\n F[(UserDemand #13, Terminal #10)] ≥ 0\n F[(Basin #12, UserDemand #13)] ≥ 0\n F[(UserDemand #6, Basin #5)] ≥ 0\n F[(LinearResistance #4, Basin #5)] ≥ 0\n F[(Basin #5, LinearResistance #4)] ≥ 0\n F[(FractionalFlow #8, Terminal #10)] ≥ 0\n F[(TabulatedRatingCurve #7, FractionalFlow #8)] ≥ 0\n F[(FractionalFlow #9, Basin #12)] ≥ 0\n F[(UserDemand #3, Basin #2)] ≥ 0\n F[(Basin #5, TabulatedRatingCurve #7)] ≥ 0\n F[(TabulatedRatingCurve #7, FractionalFlow #9)] ≥ 0\n F[(Basin #2, UserDemand #3)] ≥ 0\n F_flow_buffer_in[TabulatedRatingCurve #7] ≥ 0\n F_flow_buffer_out[TabulatedRatingCurve #7] ≥ 0",
     "crumbs": [
       "Concepts",
       "Implementation",
@@ -211,7 +211,7 @@
     "href": "concept/modelconcept.html#water-balance-equations",
     "title": "Model concept",
     "section": "1.1 Water balance equations",
-    "text": "1.1 Water balance equations\nThe water balance equation for a drainage basin (Wikipedia contributors 2022) can be defined by a first-order ordinary differential equation (ODE), where the change of the storage \\(S\\) over time is determined by the inflow fluxes minus the outflow fluxes.\n\\[\n\\frac{\\mathrm{d}S}{\\mathrm{d}t} = Q_{in} - Q_{out}\n\\]\nWe can split out the fluxes into separate terms, such as precipitation \\(P\\), evapotranspiration \\(ET\\) and runoff \\(R\\). For now other fluxes are combined into \\(Q_{rest}\\). If we define all fluxes entering our reservoir as positive, and those leaving the system as negative, all fluxes can be summed up.\n\\[\n\\frac{\\mathrm{d}S}{\\mathrm{d}t} = R + P + ET + Q_{rest}\n\\]",
+    "text": "1.1 Water balance equations\nThe water balance equation for a drainage basin (Wikipedia contributors 2022) can be defined by a first-order ordinary differential equation (ODE), where the change of the storage \\(S\\) over time is determined by the inflow fluxes minus the outflow fluxes.\n\\[\n\\frac{\\mathrm{d}S}{\\mathrm{d}t} = Q_{in} - Q_{out}\n\\]\nWe can split out the fluxes into separate terms, such as precipitation \\(P\\) and evapotranspiration \\(ET\\). For now other fluxes are combined into \\(Q_{rest}\\). If we define all fluxes entering our reservoir as positive, and those leaving the system as negative, all fluxes can be summed up.\n\\[\n\\frac{\\mathrm{d}S}{\\mathrm{d}t} = P + ET + Q_{rest}\n\\]",
     "crumbs": [
       "Concepts",
       "Numerics",
@@ -527,7 +527,7 @@
     "href": "reference/node/basin.html#static",
     "title": "Basin",
     "section": "1.1 Static",
-    "text": "1.1 Static\nThe Basin / static table can be used to set the static value of variables. The time table has a similar schema, with the time column added. A static value for a variable is only used if there is no dynamic forcing data for that variable. Specifically, if there is either no time table, it is empty, or all timestamps of that variable are missing.\n\n\n\ncolumn\ntype\nunit\nrestriction\n\n\n\n\nnode_id\nInt32\n-\nsorted\n\n\nprecipitation\nFloat64\n\\(m s^{-1}\\)\nnon-negative\n\n\npotential_evaporation\nFloat64\n\\(m s^{-1}\\)\nnon-negative\n\n\ndrainage\nFloat64\n\\(m^3 s^{-1}\\)\nnon-negative\n\n\ninfiltration\nFloat64\n\\(m^3 s^{-1}\\)\nnon-negative\n\n\nurban_runoff\nFloat64\n\\(m^3 s^{-1}\\)\nnon-negative\n\n\n\nNote that if variables are not set in the static table, default values are used when possible. These are generally zero, e.g. no precipitation, no inflow. If it is not possible to have a reasonable and safe default, a value must be provided in the static table.",
+    "text": "1.1 Static\nThe Basin / static table can be used to set the static value of variables. The time table has a similar schema, with the time column added. A static value for a variable is only used if there is no dynamic forcing data for that variable. Specifically, if there is either no time table, it is empty, or all timestamps of that variable are missing.\n\n\n\ncolumn\ntype\nunit\nrestriction\n\n\n\n\nnode_id\nInt32\n-\nsorted\n\n\nprecipitation\nFloat64\n\\(m s^{-1}\\)\nnon-negative\n\n\npotential_evaporation\nFloat64\n\\(m s^{-1}\\)\nnon-negative\n\n\ndrainage\nFloat64\n\\(m^3 s^{-1}\\)\nnon-negative\n\n\ninfiltration\nFloat64\n\\(m^3 s^{-1}\\)\nnon-negative\n\n\n\nNote that if variables are not set in the static table, default values are used when possible. These are generally zero, e.g. no precipitation, no inflow. If it is not possible to have a reasonable and safe default, a value must be provided in the static table.",
     "crumbs": [
       "Reference",
       "Nodes",