Ribasim.config is a submodule of Ribasim to handle the configuration of a Ribasim model. It is implemented using the Configurations package. A full configuration is represented by Config, which is the main API. Ribasim.config is a submodule mainly to avoid name clashes between the configuration sections and the rest of Ribasim.
Object for all information about allocation allocationnetworkids: The unique sorted allocation network IDs allocation models: The allocation models for the main network and subnetworks corresponding to allocationnetworkids mainnetworkconnections: (fromid, toid) from the main network to the subnetwork per subnetwork subnetwork_demands: The demand of an edge from the main network to a subnetwork record: A record of all flows computed by allocation optimization, eventually saved to output file
Store information for a subnetwork used for allocation.
objectivetype: The name of the type of objective used allocationnetworkid: The ID of this allocation network capacity: The capacity per edge of the allocation graph, as constrained by nodes that have a maxflowrate problem: The JuMP.jl model for solving the allocation problem Δtallocation: The time interval between consecutive allocation solves
config: The model configuration with allocation configuration in config.allocation p: Ribasim problem parameters Δt_allocation: The timestep between successive allocation solves
nodeid: node ID of the DiscreteControl node; these are not unique but repeated by the amount of conditions of this DiscreteControl node listennodeid: the ID of the node being condition on variable: the name of the variable in the condition greaterthan: The threshold value in the condition conditionvalue: The current value of each condition controlstate: Dictionary: node ID => (control state, control state start) logic_mapping: Dictionary: (control node ID, truth state) => control state record: Namedtuple with discrete control information for results
Type for storing metadata of edges in the graph: id: ID of the edge (only used for labeling flow output) type: type of the edge allocationnetworkidsource: ID of allocation network where this edge is a source (0 if not a source) fromid: the node ID of the source node toid: the node ID of the destination node allocationflow: whether this edge has a flow in an allocation graph node_ids: if this edge has allocation flow, these are all the nodes from the physical layer this edge consists of
The Model struct is an initialized model, combined with the Config used to create it and saved results. The Basic Model Interface (BMI) is implemented on the Model. A Model can be created from the path to a TOML configuration file, or a Config object.
nodeid: node ID of the Outlet node active: whether this node is active and thus contributes flow flowrate: target flow rate minflowrate: The minimal flow rate of the outlet maxflowrate: The maximum flow rate of the outlet controlmapping: dictionary from (nodeid, controlstate) to target flow rate ispid_controlled: whether the flow rate of this outlet is governed by PID control
PID control currently only supports regulating basin levels.
nodeid: node ID of the PidControl node active: whether this node is active and thus sets flow rates listennodeid: the id of the basin being controlled pidparams: a vector interpolation for parameters changing over time. The parameters are respectively target, proportional, integral, derivative, where the last three are the coefficients for the PID equation. error: the current error; basintarget - currentlevel
nodeid: node ID of the Pump node active: whether this node is active and thus contributes flow flowrate: target flow rate minflowrate: The minimal flow rate of the pump maxflowrate: The maximum flow rate of the pump controlmapping: dictionary from (nodeid, controlstate) to target flow rate ispid_controlled: whether the flow rate of this pump is governed by PID control
Rating curve from level to flow rate. The rating curve is a lookup table with linear interpolation in between. Relation can be updated in time, which is done by moving data from the time field into the tables, which is done in the update_tabulated_rating_curve callback.
Type parameter C indicates the content backing the StructVector, which can be a NamedTuple of Vectors or Arrow Primitives, and is added to avoid type instabilities.
nodeid: node ID of the TabulatedRatingCurve node active: whether this node is active and thus contributes flows tables: The current Q(h) relationships time: The time table used for updating the tables controlmapping: dictionary from (nodeid, controlstate) to Q(h) and/or active state
demand: water flux demand of user per priority over time active: whether this node is active and thus demands water allocated: water flux currently allocated to user per priority returnfactor: the factor in [0,1] of how much of the abstracted water is given back to the system minlevel: The level of the source basin below which the user does not abstract priorities: All used priority values. Each user has a demand for all these priorities, which is 0.0 if it is not provided explicitly. record: Collected data of allocation optimizations for output file.
Parse a TOML file to a Config. Keys can be overruled using keyword arguments. To overrule keys from a subsection, e.g. dt from the solver section, use underscores: solver_dt.
Add the flow capacity constraints to the allocation problem. Only finite capacities get a constraint. The constraint indices are (edgesourceid, edgedstid).
Add the source constraints to the allocation problem. The actual threshold values will be set before each allocation solve. The constraint indices are (edgesourceid, edgedstid).
Constraint: flow over source edge <= source flow in subnetwork
Certain allocation distribution types use absolute values in the objective function. Since most optimization packages do not support the absolute value function directly, New variables are introduced that act as the absolute value of an expression by posing the appropriate constraints.
Add the flow variables F to the allocation problem. The variable indices are (edgesourceid, edgedstid). Non-negativivity constraints are also immediately added to the flow variables.
Update the allocation optimization problem for the given subnetwork with the problem state and flows, solve the allocation problem and assign the results to the users.
Create the different callbacks that are used to store results and feed the simulation with new data. The different callbacks are combined to a CallbackSet that goes to the integrator. Returns the CallbackSet and the SavedValues for flow.
Convert a Real that represents the seconds passed since the simulation start to the nearest DateTime. This is used to convert between the solver’s inner float time, and the calendar.
For an element id and a vector of elements ids, get the range of indices of the last consecutive block of id. Returns the empty range 1:0 if id is not in ids.
Find the index of element x in a sorted collection a. Returns the index of x if it exists, or nothing if it doesn’t. If x occurs more than once, throw an error.
Get a sparse matrix whose sparsity matches the sparsity of the Jacobian of the ODE problem. All nodes are taken into consideration, also the ones that are inactive.
In Ribasim the Jacobian is typically sparse because each state only depends on a small number of other states.
Note: the name ‘prototype’ does not mean this code is a prototype, it comes from the naming convention of this sparsity structure in the differentialequations.jl docs.
Get the current water level of a node ID. The ID can belong to either a Basin or a LevelBoundary. storage: tells ForwardDiff whether this call is for differentiation or not
is_current_module(log::LogMessageType)::BoolReturns trueif the log message is from the current module or a submodule.See https://github.com/JuliaLogging/LoggingExtras.jl/blob/d35e7c8cfc197853ee336ace17182e6ed36dca24/src/CompositionalLoggers/earlyfiltered.jl#L39for the information available in log.
Load data from Arrow files if available, otherwise the database. Always returns a StructVector of the given struct type T, which is empty if the table is not found. This function validates the schema, and enforces the required sort order.
This is the main entry point of the application. Performs argument parsing and sets up logging for both terminal and file. Calls Ribasim.run() and handles exceptions to convert to exit codes.
Process the data in the static and time tables for a given node type. The ‘defaults’ named tuple dictates how missing data is filled in. ‘time_interpolatables’ is a vector of Symbols of parameter names for which a time interpolation (linear) object must be constructed. The control mapping for DiscreteControl is also constructed in this function. This function currently does not support node states that are defined by more than one row in a table, as is the case for TabulatedRatingCurve.
Run a Model, given a path to a TOML configuration file, or a Config object. Running a model includes initialization, solving to the end with [solve!](@ref) and writing results with write_results.
Convert a DateTime to a float that is the number of seconds since the start of the simulation. This is used to convert between the solver’s inner float time, and the calendar.
From a timeseries table time, load the most recent applicable data into table. table must be a NamedTuple of vectors with all variables that must be loaded. The most recent applicable data is non-NaN data for a given ID that is on or before t.
Update table at row index i, with the values of a given row. table must be a NamedTuple of vectors with all variables that must be loaded. The row must contain all the column names that are present in the table. If a value is missing, it is not set.
The controlled basin affects itself and the basins upstream and downstream of the controlled pump affect eachother if there is a basin upstream of the pump. The state for the integral term and the controlled basin affect eachother, and the same for the integral state and the basin upstream of the pump if it is indeed a basin.
If both the unique node upstream and the unique node downstream of these nodes are basins, then these directly depend on eachother and affect the Jacobian 2x Basins always depend on themselves.
If both the unique node upstream and the nodes down stream (or one node further if a fractional flow is in between) are basins, then the downstream basin depends on the upstream basin(s) and affect the Jacobian as many times as there are downstream basins Upstream basins always depend on themselves.
Check that nodes that have fractional flow outneighbors do not have any other type of outneighbor, that the fractions leaving a node add up to ≈1 and that the fractions are non-negative.
\[
- Q_P = P \cdot A(u).
+ Q_P = P \cdot A.
\tag{2}\]
-
Here \(P = P(t)\) is the precipitation rate and \(A\) is the wetted area. \(A\) is a function of the storage \(u = u(t)\): as the volume of water changes, the area of the free water surface may change as well, depending on the slopes of the surface waters.
+
Here \(P = P(t)\) is the precipitation rate and \(A\) is the maximum area given in the Basin / profile table. Precipitation in the Basin area is assumed to be directly added to the Basin storage. The modeler needs to ensure all precipitation enters the model, and there is no overlap in the maximum profile areas, else extra water is created. If a part of the catchment is not in any Basin profile, the modeler has to verify that water source is not forgotten. It can for instance be converted to a flow rate and added to a Basin as a FlowBoundary.
2.3 Evaporation
@@ -460,14 +460,14 @@
2.4 Infiltration and Drainage
Infiltration is provided as a lump sum for the basin. If Ribasim is coupled with MODFLOW 6, the infiltration is computed as the sum of all positive flows of the MODFLOW 6 boundary conditions in the basin:
Where \(i\) is the index of the boundary condition, \(j\) the MODFLOW 6 cell index, \(n\) the number of boundary conditions, and \(m\) the number of MODFLOW 6 cells in the basin. \(Q_{\mathrm{mf6}_{i,j}}\) is the flow computed by MODFLOW 6 for cell \(j\) for boundary condition \(i\).
Drainage is a lump sump for the basin, and consists of the sum of the absolute value of all negative flows of the MODFLOW 6 boundary conditions in the basin.
The interaction with MODFLOW 6 boundary conditions is explained in greater detail in the the MODFLOW coupling section of the documentation.
@@ -523,7 +523,7 @@
A LinearResistance connects two basins together. The flow between the two basins is determined by a linear relationship:
\[
Q = \frac{h_a - h_b}{R}
-\tag{5}\]
+\tag{6}\]
Here \(h_a\) is the water level in the first basin, \(h_b\) is the water level in the second basin, and \(R\) is the resistance of the link. A LinearResistance makes no assumptions about the direction of the flow: water flows from high to low.
@@ -671,10 +671,10 @@
3 User allocation
4 PID controller
The PID controller continuously sets the flow rate of a pump or outlet to bring the level of a certain basin closer to its setpoint. If we denote the setpoint by \(\text{SP}(t)\) and the basin level by \(y(t)\), then the error is given by \[
e(t) = \text{SP}(t) - y(t).
-\tag{6}\]
+\tag{7}\]
The output of the PID controller for the flow rate of the pump or outlet is then given by \[
Q_\text{PID}(t) = K_p e(t) + K_i\int_{t_0}^t e(\tau)\text{d}\tau + K_d \frac{\text{d}e}{\text{d}t},
-\tag{7}\]
+\tag{8}\]
for given constant parameters \(K_p,K_i,K_d\). The pump or outlet can have associated minimum and maximum flow rates \(Q_{\min}, Q_{\max}\), and so \[
Q_\text{pump/outlet} = \text{clip}(\Phi Q_\text{PID}; Q_{\min}, Q_{\max}).
\]
@@ -710,18 +710,18 @@
4 PID controller<
4.1 The derivative term
When \(K_d \ne 0\) this adds a level of complexity. We can see this by looking at the error derivative more closely: \[
\frac{\text{d}e}{\text{d}t} = \frac{\text{d}\text{SP}}{\text{d}t} - \frac{1}{A(u_\text{PID})}\frac{\text{d}u_\text{PID}}{\text{d}t},
-\] where \(A(u_\text{PID})\) is the area of the controlled basin as a function of the storage of the controlled basin \(u_\text{PID}\). The complexity arises from the fact that \(Q_\text{PID}\) is a contribution to \(\frac{\text{d}u_\text{PID}}{\text{d}t} = f_\text{PID}\), which makes Equation 7 an implicit equation for \(Q_\text{PID}\). We define
+\] where \(A(u_\text{PID})\) is the area of the controlled basin as a function of the storage of the controlled basin \(u_\text{PID}\). The complexity arises from the fact that \(Q_\text{PID}\) is a contribution to \(\frac{\text{d}u_\text{PID}}{\text{d}t} = f_\text{PID}\), which makes Equation 8 an implicit equation for \(Q_\text{PID}\). We define
that is, \(\hat{f}_\text{PID}\) is the right hand side of the ODE for the controlled basin storage state without the contribution of the PID controlled pump. The plus sign holds for an outlet and the minus sign for a pump, dictated by the way the pump and outlet connectivity to the controlled basin is enforced.
-
Using this, solving Equation 7 for \(Q_\text{PID}\) yields \[
+
Using this, solving Equation 8 for \(Q_\text{PID}\) yields \[
Q_\text{pump/outlet} = \text{clip}\left(\phi(u_\text{us})\frac{K_pe + K_iI + K_d \left(\frac{\text{d}\text{SP}}{\text{d}t}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right)}{1\pm\phi(u_\text{us})\frac{K_d}{A(u_\text{PID})}};Q_{\min},Q_{\max}\right),
\] where the clipping is again done last. Note that to compute this, \(\hat{f}_\text{PID}\) has to be known first, meaning that the PID controlled pump/outlet flow rate has to be computed after all other contributions to the PID controlled basin’s storage are known.
4.2 The sign of the parameters
-
Note by Equation 6 that the error is positive if the setpoint is larger than the basin level and negative if the setpoint is smaller than the basin level.
+
Note by Equation 7 that the error is positive if the setpoint is larger than the basin level and negative if the setpoint is smaller than the basin level.
We enforce the convention that when a pump is controlled, its edge points away from the basin, and when an outlet is controlled, its edge points towards the basin, so that the main flow direction along these edges is positive. Therefore, positive flows of the pump and outlet have opposite effects on the basin, and thus the parameters \(K_p,K_i,K_d\) of the pump and outlet must have oppositive signs to achieve the same goal.
2 Update the basi
ax = df_flow.pivot_table(index="time", columns="edge", values="flow_m3d").plot()ax.legend(bbox_to_anchor=(1.3, 1), title="Edge")
-
<matplotlib.legend.Legend at 0x7fabd80d6790>
+
<matplotlib.legend.Legend at 0x7fdc7f520e10>
diff --git a/search.json b/search.json
index d33b2fc49..95226742b 100644
--- a/search.json
+++ b/search.json
@@ -326,7 +326,7 @@
"href": "python/examples.html",
"title": "Examples",
"section": "",
- "text": "1 Basic model with static forcing\n\nfrom pathlib import Path\n\nimport geopandas as gpd\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport ribasim\n\nSetup the basins:\n\nprofile = pd.DataFrame(\n data={\n \"node_id\": [1, 1, 3, 3, 6, 6, 9, 9],\n \"area\": [0.01, 1000.0] * 4,\n \"level\": [0.0, 1.0] * 4,\n }\n)\n\n# Convert steady forcing to m/s\n# 2 mm/d precipitation, 1 mm/d evaporation\nseconds_in_day = 24 * 3600\nprecipitation = 0.002 / seconds_in_day\nevaporation = 0.001 / seconds_in_day\n\nstatic = pd.DataFrame(\n data={\n \"node_id\": [0],\n \"drainage\": [0.0],\n \"potential_evaporation\": [evaporation],\n \"infiltration\": [0.0],\n \"precipitation\": [precipitation],\n \"urban_runoff\": [0.0],\n }\n)\nstatic = static.iloc[[0, 0, 0, 0]]\nstatic[\"node_id\"] = [1, 3, 6, 9]\n\nbasin = ribasim.Basin(profile=profile, static=static)\n\nSetup linear resistance:\n\nlinear_resistance = ribasim.LinearResistance(\n static=pd.DataFrame(\n data={\"node_id\": [10, 12], \"resistance\": [5e3, (3600.0 * 24) / 100.0]}\n )\n)\n\nSetup Manning resistance:\n\nmanning_resistance = ribasim.ManningResistance(\n static=pd.DataFrame(\n data={\n \"node_id\": [2],\n \"length\": [900.0],\n \"manning_n\": [0.04],\n \"profile_width\": [6.0],\n \"profile_slope\": [3.0],\n }\n )\n)\n\nSet up a rating curve node:\n\n# Discharge: lose 1% of storage volume per day at storage = 1000.0.\nq1000 = 1000.0 * 0.01 / seconds_in_day\n\nrating_curve = ribasim.TabulatedRatingCurve(\n static=pd.DataFrame(\n data={\n \"node_id\": [4, 4],\n \"level\": [0.0, 1.0],\n \"flow_rate\": [0.0, q1000],\n }\n )\n)\n\nSetup fractional flows:\n\nfractional_flow = ribasim.FractionalFlow(\n static=pd.DataFrame(\n data={\n \"node_id\": [5, 8, 13],\n \"fraction\": [0.3, 0.6, 0.1],\n }\n )\n)\n\nSetup pump:\n\npump = ribasim.Pump(\n static=pd.DataFrame(\n data={\n \"node_id\": [7],\n \"flow_rate\": [0.5 / 3600],\n }\n )\n)\n\nSetup level boundary:\n\nlevel_boundary = ribasim.LevelBoundary(\n static=pd.DataFrame(\n data={\n \"node_id\": [11, 17],\n \"level\": [0.5, 1.5],\n }\n )\n)\n\nSetup flow boundary:\n\nflow_boundary = ribasim.FlowBoundary(\n static=pd.DataFrame(\n data={\n \"node_id\": [15, 16],\n \"flow_rate\": [1e-4, 1e-4],\n }\n )\n)\n\nSetup terminal:\n\nterminal = ribasim.Terminal(\n static=pd.DataFrame(\n data={\n \"node_id\": [14],\n }\n )\n)\n\nSet up the nodes:\n\nxy = np.array(\n [\n (0.0, 0.0), # 1: Basin,\n (1.0, 0.0), # 2: ManningResistance\n (2.0, 0.0), # 3: Basin\n (3.0, 0.0), # 4: TabulatedRatingCurve\n (3.0, 1.0), # 5: FractionalFlow\n (3.0, 2.0), # 6: Basin\n (4.0, 1.0), # 7: Pump\n (4.0, 0.0), # 8: FractionalFlow\n (5.0, 0.0), # 9: Basin\n (6.0, 0.0), # 10: LinearResistance\n (2.0, 2.0), # 11: LevelBoundary\n (2.0, 1.0), # 12: LinearResistance\n (3.0, -1.0), # 13: FractionalFlow\n (3.0, -2.0), # 14: Terminal\n (3.0, 3.0), # 15: FlowBoundary\n (0.0, 1.0), # 16: FlowBoundary\n (6.0, 1.0), # 17: LevelBoundary\n ]\n)\nnode_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n\nnode_id, node_type = ribasim.Node.node_ids_and_types(\n basin,\n manning_resistance,\n rating_curve,\n pump,\n fractional_flow,\n linear_resistance,\n level_boundary,\n flow_boundary,\n terminal,\n)\n\n# Make sure the feature id starts at 1: explicitly give an index.\nnode = ribasim.Node(\n df=gpd.GeoDataFrame(\n data={\"type\": node_type},\n index=pd.Index(node_id, name=\"fid\"),\n geometry=node_xy,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the edges:\n\nfrom_id = np.array(\n [1, 2, 3, 4, 4, 5, 6, 8, 7, 9, 11, 12, 4, 13, 15, 16, 10], dtype=np.int64\n)\nto_id = np.array(\n [2, 3, 4, 5, 8, 6, 7, 9, 9, 10, 12, 3, 13, 14, 6, 1, 17], dtype=np.int64\n)\nlines = node.geometry_from_connectivity(from_id, to_id)\nedge = ribasim.Edge(\n df=gpd.GeoDataFrame(\n data={\n \"from_node_id\": from_id,\n \"to_node_id\": to_id,\n \"edge_type\": len(from_id) * [\"flow\"],\n },\n geometry=lines,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup a model:\n\nmodel = ribasim.Model(\n network=ribasim.Network(\n node=node,\n edge=edge,\n ),\n basin=basin,\n level_boundary=level_boundary,\n flow_boundary=flow_boundary,\n pump=pump,\n linear_resistance=linear_resistance,\n manning_resistance=manning_resistance,\n tabulated_rating_curve=rating_curve,\n fractional_flow=fractional_flow,\n terminal=terminal,\n starttime=\"2020-01-01 00:00:00\",\n endtime=\"2021-01-01 00:00:00\",\n)\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n<Axes: >\n\n\n\n\n\nWrite the model to a TOML and GeoPackage:\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"basic/ribasim.toml\")\n\nPosixPath('data/basic/ribasim.toml')\n\n\n\n\n2 Update the basic model with transient forcing\nThis assumes you have already created the basic model with static forcing.\n\nimport numpy as np\nimport pandas as pd\nimport ribasim\nimport xarray as xr\n\n\nmodel = ribasim.Model(filepath=datadir / \"basic/ribasim.toml\")\n\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\nWe’ll use xarray to easily broadcast the values.\n\ntimeseries = (\n pd.DataFrame(\n data={\n \"node_id\": 1,\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 )\n .set_index(\"time\")\n .to_xarray()\n)\n\nbasin_ids = model.basin.static.df[\"node_id\"].to_numpy()\nbasin_nodes = xr.DataArray(\n np.ones(len(basin_ids)), coords={\"node_id\": basin_ids}, dims=[\"node_id\"]\n)\nforcing = (timeseries * basin_nodes).to_dataframe().reset_index()\n\n\nstate = pd.DataFrame(\n data={\n \"node_id\": basin_ids,\n \"level\": 1.4,\n }\n)\n\n\nmodel.basin.time.df = forcing\nmodel.basin.state.df = state\n\n\nmodel.write(datadir / \"basic_transient/ribasim.toml\")\n\nPosixPath('data/basic_transient/ribasim.toml')\n\n\nNow run the model with ribasim basic_transient/ribasim.toml. After running the model, read back the results:\n\ndf_basin = pd.read_feather(datadir / \"basic_transient/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<Axes: xlabel='time'>\n\n\n\n\n\n\ndf_flow = pd.read_feather(datadir / \"basic_transient/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 * 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<matplotlib.legend.Legend at 0x7fabd80d6790>\n\n\n\n\n\n\ntype(df_flow)\n\npandas.core.frame.DataFrame\n\n\n\n\n3 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 around a target level (setpoint) by two connected pumps. These two pumps 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.\nSet up the nodes:\n\nxy = np.array(\n [\n (0.0, 0.0), # 1: Basin\n (1.0, 1.0), # 2: Pump\n (1.0, -1.0), # 3: Pump\n (2.0, 0.0), # 4: LevelBoundary\n (-1.0, 0.0), # 5: TabulatedRatingCurve\n (-2.0, 0.0), # 6: Terminal\n (1.0, 0.0), # 7: DiscreteControl\n ]\n)\n\nnode_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n\nnode_type = [\n \"Basin\",\n \"Pump\",\n \"Pump\",\n \"LevelBoundary\",\n \"TabulatedRatingCurve\",\n \"Terminal\",\n \"DiscreteControl\",\n]\n\n# Make sure the feature id starts at 1: explicitly give an index.\nnode = ribasim.Node(\n df=gpd.GeoDataFrame(\n data={\"type\": node_type},\n index=pd.Index(np.arange(len(xy)) + 1, name=\"fid\"),\n geometry=node_xy,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the edges:\n\nfrom_id = np.array([1, 3, 4, 2, 1, 5, 7, 7], dtype=np.int64)\nto_id = np.array([3, 4, 2, 1, 5, 6, 2, 3], dtype=np.int64)\n\nedge_type = 6 * [\"flow\"] + 2 * [\"control\"]\n\nlines = node.geometry_from_connectivity(from_id, to_id)\nedge = ribasim.Edge(\n df=gpd.GeoDataFrame(\n data={\"from_node_id\": from_id, \"to_node_id\": to_id, \"edge_type\": edge_type},\n geometry=lines,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the basins:\n\nprofile = pd.DataFrame(\n data={\n \"node_id\": [1, 1],\n \"area\": [1000.0, 1000.0],\n \"level\": [0.0, 1.0],\n }\n)\n\nstatic = pd.DataFrame(\n data={\n \"node_id\": [1],\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\nstate = pd.DataFrame(data={\"node_id\": [1], \"level\": [20.0]})\n\nbasin = ribasim.Basin(profile=profile, static=static, state=state)\n\nSetup the discrete control:\n\ncondition = pd.DataFrame(\n data={\n \"node_id\": 3 * [7],\n \"listen_feature_id\": 3 * [1],\n \"variable\": 3 * [\"level\"],\n \"greater_than\": [5.0, 10.0, 15.0], # min, setpoint, max\n }\n)\n\nlogic = pd.DataFrame(\n data={\n \"node_id\": 5 * [7],\n \"truth_state\": [\"FFF\", \"U**\", \"T*F\", \"**D\", \"TTT\"],\n \"control_state\": [\"in\", \"in\", \"none\", \"out\", \"out\"],\n }\n)\n\ndiscrete_control = ribasim.DiscreteControl(condition=condition, logic=logic)\n\nThe above control logic can be summarized as follows: - If the level gets above the maximum, activate the control state “out” until the setpoint is reached; - If the level gets below the minimum, active the control state “in” until the setpoint is reached; - Otherwise activate the control state “none”.\nSetup the pump:\n\npump = ribasim.Pump(\n static=pd.DataFrame(\n data={\n \"node_id\": 3 * [2] + 3 * [3],\n \"control_state\": 2 * [\"none\", \"in\", \"out\"],\n \"flow_rate\": [0.0, 2e-3, 0.0, 0.0, 0.0, 2e-3],\n }\n )\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\nlevel_boundary = ribasim.LevelBoundary(\n static=pd.DataFrame(data={\"node_id\": [4], \"level\": [10.0]})\n)\n\nSetup the rating curve:\n\nrating_curve = ribasim.TabulatedRatingCurve(\n static=pd.DataFrame(\n data={\"node_id\": 2 * [5], \"level\": [2.0, 15.0], \"flow_rate\": [0.0, 1e-3]}\n )\n)\n\nSetup the terminal:\n\nterminal = ribasim.Terminal(static=pd.DataFrame(data={\"node_id\": [6]}))\n\nSetup a model:\n\nmodel = ribasim.Model(\n network=ribasim.Network(\n node=node,\n edge=edge,\n ),\n basin=basin,\n pump=pump,\n level_boundary=level_boundary,\n tabulated_rating_curve=rating_curve,\n terminal=terminal,\n discrete_control=discrete_control,\n starttime=\"2020-01-01 00:00:00\",\n endtime=\"2021-01-01 00:00:00\",\n)\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n<Axes: >\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_setpoint_with_minmax/ribasim.toml\")\n\nPosixPath('data/level_setpoint_with_minmax/ribasim.toml')\n\n\nNow run the model with level_setpoint_with_minmax/ribasim.toml. After running the model, read back the results:\n\nfrom matplotlib.dates import date2num\n\ndf_basin = pd.read_feather(datadir / \"level_setpoint_with_minmax/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\ndf_control = pd.read_feather(\n datadir / \"level_setpoint_with_minmax/results/control.arrow\"\n)\n\ny_min, y_max = ax.get_ybound()\nax.fill_between(df_control.time[:2], 2 * [y_min], 2 * [y_max], alpha=0.2, color=\"C0\")\nax.fill_between(df_control.time[2:4], 2 * [y_min], 2 * [y_max], alpha=0.2, color=\"C0\")\n\nax.set_xticks(\n date2num(df_control.time).tolist(),\n df_control.control_state.tolist(),\n rotation=50,\n)\n\nax.set_yticks(greater_than, [\"min\", \"setpoint\", \"max\"])\nax.set_ylabel(\"level\")\nplt.show()\n\n\n\n\nThe highlighted regions show where a pump is active.\nLet’s print an overview of what happened with control:\n\nmodel.print_discrete_control_record(\n datadir / \"level_setpoint_with_minmax/results/control.arrow\"\n)\n\n0. At 2020-01-01 00:00:00 the control node with ID 7 reached truth state TTT:\n For node ID 1 (Basin): level > 5.0\n For node ID 1 (Basin): level > 10.0\n For node ID 1 (Basin): level > 15.0\n\n This yielded control state \"out\":\n For node ID 2 (Pump): active = <NA>, flow_rate = 0.0, min_flow_rate = nan, max_flow_rate = nan\n For node ID 3 (Pump): active = <NA>, flow_rate = 0.002, min_flow_rate = nan, max_flow_rate = nan\n\n1. At 2020-02-08 19:02:21.861000 the control node with ID 7 reached truth state TFF:\n For node ID 1 (Basin): level > 5.0\n For node ID 1 (Basin): level < 10.0\n For node ID 1 (Basin): level < 15.0\n\n This yielded control state \"none\":\n For node ID 2 (Pump): active = <NA>, flow_rate = 0.0, min_flow_rate = nan, max_flow_rate = nan\n For node ID 3 (Pump): active = <NA>, flow_rate = 0.0, min_flow_rate = nan, max_flow_rate = nan\n\n2. At 2020-07-05 08:56:10.319000 the control node with ID 7 reached truth state FFF:\n For node ID 1 (Basin): level < 5.0\n For node ID 1 (Basin): level < 10.0\n For node ID 1 (Basin): level < 15.0\n\n This yielded control state \"in\":\n For node ID 2 (Pump): active = <NA>, flow_rate = 0.002, min_flow_rate = nan, max_flow_rate = nan\n For node ID 3 (Pump): active = <NA>, flow_rate = 0.0, min_flow_rate = nan, max_flow_rate = nan\n\n3. At 2020-08-11 06:05:15.592000 the control node with ID 7 reached truth state TTF:\n For node ID 1 (Basin): level > 5.0\n For node ID 1 (Basin): level > 10.0\n For node ID 1 (Basin): level < 15.0\n\n This yielded control state \"none\":\n For node ID 2 (Pump): active = <NA>, flow_rate = 0.0, min_flow_rate = nan, max_flow_rate = nan\n For node ID 3 (Pump): active = <NA>, flow_rate = 0.0, min_flow_rate = nan, max_flow_rate = nan\n\n\n\nNote that crossing direction specific truth states (containing “U”, “D”) are not present in this overview even though they are part of the control logic. This is because in the control logic for this model these truth states are only used to sustain control states, while the overview only shows changes in control states.\n\n\n4 Model with PID control\nSet up the nodes:\n\nxy = np.array(\n [\n (0.0, 0.0), # 1: FlowBoundary\n (1.0, 0.0), # 2: Basin\n (2.0, 0.5), # 3: Pump\n (3.0, 0.0), # 4: LevelBoundary\n (1.5, 1.0), # 5: PidControl\n (2.0, -0.5), # 6: outlet\n (1.5, -1.0), # 7: PidControl\n ]\n)\n\nnode_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n\nnode_type = [\n \"FlowBoundary\",\n \"Basin\",\n \"Pump\",\n \"LevelBoundary\",\n \"PidControl\",\n \"Outlet\",\n \"PidControl\",\n]\n\n# Make sure the feature id starts at 1: explicitly give an index.\nnode = ribasim.Node(\n df=gpd.GeoDataFrame(\n data={\"type\": node_type},\n index=pd.Index(np.arange(len(xy)) + 1, name=\"fid\"),\n geometry=node_xy,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the edges:\n\nfrom_id = np.array([1, 2, 3, 4, 6, 5, 7], dtype=np.int64)\nto_id = np.array([2, 3, 4, 6, 2, 3, 6], dtype=np.int64)\n\nlines = node.geometry_from_connectivity(from_id, to_id)\nedge = ribasim.Edge(\n df=gpd.GeoDataFrame(\n data={\n \"from_node_id\": from_id,\n \"to_node_id\": to_id,\n \"edge_type\": 5 * [\"flow\"] + 2 * [\"control\"],\n },\n geometry=lines,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the basins:\n\nprofile = pd.DataFrame(\n data={\"node_id\": [2, 2], \"level\": [0.0, 1.0], \"area\": [1000.0, 1000.0]}\n)\n\nstatic = pd.DataFrame(\n data={\n \"node_id\": [2],\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\nstate = pd.DataFrame(\n data={\n \"node_id\": [2],\n \"level\": [6.0],\n }\n)\n\nbasin = ribasim.Basin(profile=profile, static=static, state=state)\n\nSetup the pump:\n\npump = ribasim.Pump(\n static=pd.DataFrame(\n data={\n \"node_id\": [3],\n \"flow_rate\": [0.0], # Will be overwritten by PID controller\n }\n )\n)\n\nSetup the outlet:\n\noutlet = ribasim.Outlet(\n static=pd.DataFrame(\n data={\n \"node_id\": [6],\n \"flow_rate\": [0.0], # Will be overwritten by PID controller\n }\n )\n)\n\nSetup flow boundary:\n\nflow_boundary = ribasim.FlowBoundary(\n static=pd.DataFrame(data={\"node_id\": [1], \"flow_rate\": [1e-3]})\n)\n\nSetup flow boundary:\n\nlevel_boundary = ribasim.LevelBoundary(\n static=pd.DataFrame(\n data={\n \"node_id\": [4],\n \"level\": [1.0], # Not relevant\n }\n )\n)\n\nSetup PID control:\n\npid_control = ribasim.PidControl(\n time=pd.DataFrame(\n data={\n \"node_id\": 4 * [5, 7],\n \"time\": [\n \"2020-01-01 00:00:00\",\n \"2020-01-01 00:00:00\",\n \"2020-05-01 00:00:00\",\n \"2020-05-01 00:00:00\",\n \"2020-07-01 00:00:00\",\n \"2020-07-01 00:00:00\",\n \"2020-12-01 00:00:00\",\n \"2020-12-01 00:00:00\",\n ],\n \"listen_node_id\": 4 * [2, 2],\n \"target\": [5.0, 5.0, 5.0, 5.0, 7.5, 7.5, 7.5, 7.5],\n \"proportional\": 4 * [-1e-3, 1e-3],\n \"integral\": 4 * [-1e-7, 1e-7],\n \"derivative\": 4 * [0.0, 0.0],\n }\n )\n)\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 a model:\n\nmodel = ribasim.Model(\n network=ribasim.Network(\n node=node,\n edge=edge,\n ),\n basin=basin,\n flow_boundary=flow_boundary,\n level_boundary=level_boundary,\n pump=pump,\n outlet=outlet,\n pid_control=pid_control,\n starttime=\"2020-01-01 00:00:00\",\n endtime=\"2020-12-01 00:00:00\",\n)\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n<Axes: >\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\ntarget_levels = model.pid_control.time.df.target.to_numpy()[:4]\ntimes = date2num(model.pid_control.time.df.time)[:4]\nax.plot(times, target_levels, color=\"k\", ls=\":\", label=\"target level\")\npass\n\n\n\n\n\n\n5 Model with allocation\nSetup the nodes:\n\nxy = np.array(\n [\n (0.0, 0.0), # 1: FlowBoundary\n (1.0, 0.0), # 2: Basin\n (1.0, 1.0), # 3: User\n (2.0, 0.0), # 4: LinearResistance\n (3.0, 0.0), # 5: Basin\n (3.0, 1.0), # 6: User\n (4.0, 0.0), # 7: TabulatedRatingCurve\n (4.5, 0.0), # 8: FractionalFlow\n (4.5, 0.5), # 9: FractionalFlow\n (5.0, 0.0), # 10: Terminal\n (4.5, 0.25), # 11: DiscreteControl\n (4.5, 1.0), # 12: Basin\n (5.0, 1.0), # 13: User\n ]\n)\nnode_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n\nnode_type = [\n \"FlowBoundary\",\n \"Basin\",\n \"User\",\n \"LinearResistance\",\n \"Basin\",\n \"User\",\n \"TabulatedRatingCurve\",\n \"FractionalFlow\",\n \"FractionalFlow\",\n \"Terminal\",\n \"DiscreteControl\",\n \"Basin\",\n \"User\",\n]\n\n# All nodes belong to allocation network id 1\nnode = ribasim.Node(\n df=gpd.GeoDataFrame(\n data={\"type\": node_type, \"allocation_network_id\": 1},\n index=pd.Index(np.arange(len(xy)) + 1, name=\"fid\"),\n geometry=node_xy,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the edges:\n\nfrom_id = np.array(\n [1, 2, 2, 4, 5, 5, 7, 3, 6, 7, 8, 9, 12, 13, 11, 11],\n dtype=np.int64,\n)\nto_id = np.array(\n [2, 3, 4, 5, 6, 7, 8, 2, 5, 9, 10, 12, 13, 10, 8, 9],\n dtype=np.int64,\n)\n# Denote the first edge, 1 => 2, as a source edge for\n# allocation network 1\nallocation_network_id = len(from_id) * [None]\nallocation_network_id[0] = 1\nlines = node.geometry_from_connectivity(from_id, to_id)\nedge = ribasim.Edge(\n df=gpd.GeoDataFrame(\n data={\n \"from_node_id\": from_id,\n \"to_node_id\": to_id,\n \"edge_type\": (len(from_id) - 2) * [\"flow\"] + 2 * [\"control\"],\n \"allocation_network_id\": allocation_network_id,\n },\n geometry=lines,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the basins:\n\nprofile = pd.DataFrame(\n data={\n \"node_id\": [2, 2, 5, 5, 12, 12],\n \"area\": 300_000.0,\n \"level\": 3 * [0.0, 1.0],\n }\n)\n\nstatic = pd.DataFrame(\n data={\n \"node_id\": [2, 5, 12],\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\nstate = pd.DataFrame(data={\"node_id\": [2, 5, 12], \"level\": 1.0})\n\nbasin = ribasim.Basin(profile=profile, static=static, state=state)\n\nSetup the flow boundary:\n\nflow_boundary = ribasim.FlowBoundary(\n static=pd.DataFrame(\n data={\n \"node_id\": [1],\n \"flow_rate\": 2.0,\n }\n )\n)\n\nSetup the linear resistance:\n\nlinear_resistance = ribasim.LinearResistance(\n static=pd.DataFrame(\n data={\n \"node_id\": [4],\n \"resistance\": 0.06,\n }\n )\n)\n\nSetup the tabulated rating curve:\n\ntabulated_rating_curve = ribasim.TabulatedRatingCurve(\n static=pd.DataFrame(\n data={\n \"node_id\": 7,\n \"level\": [0.0, 0.5, 1.0],\n \"flow_rate\": [0.0, 0.0, 2.0],\n }\n )\n)\n\nSetup the fractional flow:\n\nfractional_flow = ribasim.FractionalFlow(\n static=pd.DataFrame(\n data={\n \"node_id\": [8, 8, 9, 9],\n \"fraction\": [0.6, 0.9, 0.4, 0.1],\n \"control_state\": [\"divert\", \"close\", \"divert\", \"close\"],\n }\n )\n)\n\nSetup the terminal:\n\nterminal = ribasim.Terminal(\n static=pd.DataFrame(\n data={\n \"node_id\": [10],\n }\n )\n)\n\nSetup the discrete control:\n\ncondition = pd.DataFrame(\n data={\n \"node_id\": [11],\n \"listen_feature_id\": 5,\n \"variable\": \"level\",\n \"greater_than\": 0.52,\n }\n)\n\nlogic = pd.DataFrame(\n data={\n \"node_id\": 11,\n \"truth_state\": [\"T\", \"F\"],\n \"control_state\": [\"divert\", \"close\"],\n }\n)\n\ndiscrete_control = ribasim.DiscreteControl(condition=condition, logic=logic)\n\nSetup the users:\n\nuser = ribasim.User(\n static=pd.DataFrame(\n data={\n \"node_id\": [6, 13],\n \"demand\": [1.5, 1.0],\n \"return_factor\": 0.0,\n \"min_level\": -1.0,\n \"priority\": [1, 3],\n }\n ),\n time=pd.DataFrame(\n data={\n \"node_id\": [3, 3, 3, 3],\n \"demand\": [0.0, 1.0, 1.2, 1.2],\n \"priority\": [1, 1, 2, 2],\n \"return_factor\": 0.0,\n \"min_level\": -1.0,\n \"time\": 2 * [\"2020-01-01 00:00:00\", \"2020-01-20 00:00:00\"],\n }\n ),\n)\n\nSetup the allocation:\n\nallocation = ribasim.Allocation(use_allocation=True, timestep=86400)\n\nSetup a model:\n\nmodel = ribasim.Model(\n network=ribasim.Network(\n node=node,\n edge=edge,\n ),\n basin=basin,\n flow_boundary=flow_boundary,\n linear_resistance=linear_resistance,\n tabulated_rating_curve=tabulated_rating_curve,\n terminal=terminal,\n user=user,\n discrete_control=discrete_control,\n fractional_flow=fractional_flow,\n allocation=allocation,\n starttime=\"2020-01-01 00:00:00\",\n endtime=\"2020-01-20 00:00:00\",\n)\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n<Axes: >\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=[\"user_node_id\", \"priority\"],\n values=[\"demand\", \"allocated\", \"abstracted\"],\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))\nfig.suptitle(f\"Objective type: {model.allocation.objective_type}\", fontweight=\"bold\")\n\ndf_allocation_wide[\"demand\"].plot(ax=axs[0], ls=\":\")\ndf_allocation_wide[\"allocated\"].plot(ax=axs[1], ls=\"--\")\ndf_allocation_wide[\"abstracted\"].plot(ax=axs[2])\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\nSome things to note about this plot:\n\nAbstraction behaves somewhat erratically at the start of the simulation. This is because allocation is based on flows computed in the physical layer, and at the start of the simulation these are not known yet.\nAlthough there is a plotted line for abstraction per priority, abstraction is actually 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]')"
+ "text": "1 Basic model with static forcing\n\nfrom pathlib import Path\n\nimport geopandas as gpd\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport ribasim\n\nSetup the basins:\n\nprofile = pd.DataFrame(\n data={\n \"node_id\": [1, 1, 3, 3, 6, 6, 9, 9],\n \"area\": [0.01, 1000.0] * 4,\n \"level\": [0.0, 1.0] * 4,\n }\n)\n\n# Convert steady forcing to m/s\n# 2 mm/d precipitation, 1 mm/d evaporation\nseconds_in_day = 24 * 3600\nprecipitation = 0.002 / seconds_in_day\nevaporation = 0.001 / seconds_in_day\n\nstatic = pd.DataFrame(\n data={\n \"node_id\": [0],\n \"drainage\": [0.0],\n \"potential_evaporation\": [evaporation],\n \"infiltration\": [0.0],\n \"precipitation\": [precipitation],\n \"urban_runoff\": [0.0],\n }\n)\nstatic = static.iloc[[0, 0, 0, 0]]\nstatic[\"node_id\"] = [1, 3, 6, 9]\n\nbasin = ribasim.Basin(profile=profile, static=static)\n\nSetup linear resistance:\n\nlinear_resistance = ribasim.LinearResistance(\n static=pd.DataFrame(\n data={\"node_id\": [10, 12], \"resistance\": [5e3, (3600.0 * 24) / 100.0]}\n )\n)\n\nSetup Manning resistance:\n\nmanning_resistance = ribasim.ManningResistance(\n static=pd.DataFrame(\n data={\n \"node_id\": [2],\n \"length\": [900.0],\n \"manning_n\": [0.04],\n \"profile_width\": [6.0],\n \"profile_slope\": [3.0],\n }\n )\n)\n\nSet up a rating curve node:\n\n# Discharge: lose 1% of storage volume per day at storage = 1000.0.\nq1000 = 1000.0 * 0.01 / seconds_in_day\n\nrating_curve = ribasim.TabulatedRatingCurve(\n static=pd.DataFrame(\n data={\n \"node_id\": [4, 4],\n \"level\": [0.0, 1.0],\n \"flow_rate\": [0.0, q1000],\n }\n )\n)\n\nSetup fractional flows:\n\nfractional_flow = ribasim.FractionalFlow(\n static=pd.DataFrame(\n data={\n \"node_id\": [5, 8, 13],\n \"fraction\": [0.3, 0.6, 0.1],\n }\n )\n)\n\nSetup pump:\n\npump = ribasim.Pump(\n static=pd.DataFrame(\n data={\n \"node_id\": [7],\n \"flow_rate\": [0.5 / 3600],\n }\n )\n)\n\nSetup level boundary:\n\nlevel_boundary = ribasim.LevelBoundary(\n static=pd.DataFrame(\n data={\n \"node_id\": [11, 17],\n \"level\": [0.5, 1.5],\n }\n )\n)\n\nSetup flow boundary:\n\nflow_boundary = ribasim.FlowBoundary(\n static=pd.DataFrame(\n data={\n \"node_id\": [15, 16],\n \"flow_rate\": [1e-4, 1e-4],\n }\n )\n)\n\nSetup terminal:\n\nterminal = ribasim.Terminal(\n static=pd.DataFrame(\n data={\n \"node_id\": [14],\n }\n )\n)\n\nSet up the nodes:\n\nxy = np.array(\n [\n (0.0, 0.0), # 1: Basin,\n (1.0, 0.0), # 2: ManningResistance\n (2.0, 0.0), # 3: Basin\n (3.0, 0.0), # 4: TabulatedRatingCurve\n (3.0, 1.0), # 5: FractionalFlow\n (3.0, 2.0), # 6: Basin\n (4.0, 1.0), # 7: Pump\n (4.0, 0.0), # 8: FractionalFlow\n (5.0, 0.0), # 9: Basin\n (6.0, 0.0), # 10: LinearResistance\n (2.0, 2.0), # 11: LevelBoundary\n (2.0, 1.0), # 12: LinearResistance\n (3.0, -1.0), # 13: FractionalFlow\n (3.0, -2.0), # 14: Terminal\n (3.0, 3.0), # 15: FlowBoundary\n (0.0, 1.0), # 16: FlowBoundary\n (6.0, 1.0), # 17: LevelBoundary\n ]\n)\nnode_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n\nnode_id, node_type = ribasim.Node.node_ids_and_types(\n basin,\n manning_resistance,\n rating_curve,\n pump,\n fractional_flow,\n linear_resistance,\n level_boundary,\n flow_boundary,\n terminal,\n)\n\n# Make sure the feature id starts at 1: explicitly give an index.\nnode = ribasim.Node(\n df=gpd.GeoDataFrame(\n data={\"type\": node_type},\n index=pd.Index(node_id, name=\"fid\"),\n geometry=node_xy,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the edges:\n\nfrom_id = np.array(\n [1, 2, 3, 4, 4, 5, 6, 8, 7, 9, 11, 12, 4, 13, 15, 16, 10], dtype=np.int64\n)\nto_id = np.array(\n [2, 3, 4, 5, 8, 6, 7, 9, 9, 10, 12, 3, 13, 14, 6, 1, 17], dtype=np.int64\n)\nlines = node.geometry_from_connectivity(from_id, to_id)\nedge = ribasim.Edge(\n df=gpd.GeoDataFrame(\n data={\n \"from_node_id\": from_id,\n \"to_node_id\": to_id,\n \"edge_type\": len(from_id) * [\"flow\"],\n },\n geometry=lines,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup a model:\n\nmodel = ribasim.Model(\n network=ribasim.Network(\n node=node,\n edge=edge,\n ),\n basin=basin,\n level_boundary=level_boundary,\n flow_boundary=flow_boundary,\n pump=pump,\n linear_resistance=linear_resistance,\n manning_resistance=manning_resistance,\n tabulated_rating_curve=rating_curve,\n fractional_flow=fractional_flow,\n terminal=terminal,\n starttime=\"2020-01-01 00:00:00\",\n endtime=\"2021-01-01 00:00:00\",\n)\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n<Axes: >\n\n\n\n\n\nWrite the model to a TOML and GeoPackage:\n\ndatadir = Path(\"data\")\nmodel.write(datadir / \"basic/ribasim.toml\")\n\nPosixPath('data/basic/ribasim.toml')\n\n\n\n\n2 Update the basic model with transient forcing\nThis assumes you have already created the basic model with static forcing.\n\nimport numpy as np\nimport pandas as pd\nimport ribasim\nimport xarray as xr\n\n\nmodel = ribasim.Model(filepath=datadir / \"basic/ribasim.toml\")\n\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\nWe’ll use xarray to easily broadcast the values.\n\ntimeseries = (\n pd.DataFrame(\n data={\n \"node_id\": 1,\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 )\n .set_index(\"time\")\n .to_xarray()\n)\n\nbasin_ids = model.basin.static.df[\"node_id\"].to_numpy()\nbasin_nodes = xr.DataArray(\n np.ones(len(basin_ids)), coords={\"node_id\": basin_ids}, dims=[\"node_id\"]\n)\nforcing = (timeseries * basin_nodes).to_dataframe().reset_index()\n\n\nstate = pd.DataFrame(\n data={\n \"node_id\": basin_ids,\n \"level\": 1.4,\n }\n)\n\n\nmodel.basin.time.df = forcing\nmodel.basin.state.df = state\n\n\nmodel.write(datadir / \"basic_transient/ribasim.toml\")\n\nPosixPath('data/basic_transient/ribasim.toml')\n\n\nNow run the model with ribasim basic_transient/ribasim.toml. After running the model, read back the results:\n\ndf_basin = pd.read_feather(datadir / \"basic_transient/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<Axes: xlabel='time'>\n\n\n\n\n\n\ndf_flow = pd.read_feather(datadir / \"basic_transient/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 * 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<matplotlib.legend.Legend at 0x7fdc7f520e10>\n\n\n\n\n\n\ntype(df_flow)\n\npandas.core.frame.DataFrame\n\n\n\n\n3 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 around a target level (setpoint) by two connected pumps. These two pumps 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.\nSet up the nodes:\n\nxy = np.array(\n [\n (0.0, 0.0), # 1: Basin\n (1.0, 1.0), # 2: Pump\n (1.0, -1.0), # 3: Pump\n (2.0, 0.0), # 4: LevelBoundary\n (-1.0, 0.0), # 5: TabulatedRatingCurve\n (-2.0, 0.0), # 6: Terminal\n (1.0, 0.0), # 7: DiscreteControl\n ]\n)\n\nnode_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n\nnode_type = [\n \"Basin\",\n \"Pump\",\n \"Pump\",\n \"LevelBoundary\",\n \"TabulatedRatingCurve\",\n \"Terminal\",\n \"DiscreteControl\",\n]\n\n# Make sure the feature id starts at 1: explicitly give an index.\nnode = ribasim.Node(\n df=gpd.GeoDataFrame(\n data={\"type\": node_type},\n index=pd.Index(np.arange(len(xy)) + 1, name=\"fid\"),\n geometry=node_xy,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the edges:\n\nfrom_id = np.array([1, 3, 4, 2, 1, 5, 7, 7], dtype=np.int64)\nto_id = np.array([3, 4, 2, 1, 5, 6, 2, 3], dtype=np.int64)\n\nedge_type = 6 * [\"flow\"] + 2 * [\"control\"]\n\nlines = node.geometry_from_connectivity(from_id, to_id)\nedge = ribasim.Edge(\n df=gpd.GeoDataFrame(\n data={\"from_node_id\": from_id, \"to_node_id\": to_id, \"edge_type\": edge_type},\n geometry=lines,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the basins:\n\nprofile = pd.DataFrame(\n data={\n \"node_id\": [1, 1],\n \"area\": [1000.0, 1000.0],\n \"level\": [0.0, 1.0],\n }\n)\n\nstatic = pd.DataFrame(\n data={\n \"node_id\": [1],\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\nstate = pd.DataFrame(data={\"node_id\": [1], \"level\": [20.0]})\n\nbasin = ribasim.Basin(profile=profile, static=static, state=state)\n\nSetup the discrete control:\n\ncondition = pd.DataFrame(\n data={\n \"node_id\": 3 * [7],\n \"listen_feature_id\": 3 * [1],\n \"variable\": 3 * [\"level\"],\n \"greater_than\": [5.0, 10.0, 15.0], # min, setpoint, max\n }\n)\n\nlogic = pd.DataFrame(\n data={\n \"node_id\": 5 * [7],\n \"truth_state\": [\"FFF\", \"U**\", \"T*F\", \"**D\", \"TTT\"],\n \"control_state\": [\"in\", \"in\", \"none\", \"out\", \"out\"],\n }\n)\n\ndiscrete_control = ribasim.DiscreteControl(condition=condition, logic=logic)\n\nThe above control logic can be summarized as follows: - If the level gets above the maximum, activate the control state “out” until the setpoint is reached; - If the level gets below the minimum, active the control state “in” until the setpoint is reached; - Otherwise activate the control state “none”.\nSetup the pump:\n\npump = ribasim.Pump(\n static=pd.DataFrame(\n data={\n \"node_id\": 3 * [2] + 3 * [3],\n \"control_state\": 2 * [\"none\", \"in\", \"out\"],\n \"flow_rate\": [0.0, 2e-3, 0.0, 0.0, 0.0, 2e-3],\n }\n )\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\nlevel_boundary = ribasim.LevelBoundary(\n static=pd.DataFrame(data={\"node_id\": [4], \"level\": [10.0]})\n)\n\nSetup the rating curve:\n\nrating_curve = ribasim.TabulatedRatingCurve(\n static=pd.DataFrame(\n data={\"node_id\": 2 * [5], \"level\": [2.0, 15.0], \"flow_rate\": [0.0, 1e-3]}\n )\n)\n\nSetup the terminal:\n\nterminal = ribasim.Terminal(static=pd.DataFrame(data={\"node_id\": [6]}))\n\nSetup a model:\n\nmodel = ribasim.Model(\n network=ribasim.Network(\n node=node,\n edge=edge,\n ),\n basin=basin,\n pump=pump,\n level_boundary=level_boundary,\n tabulated_rating_curve=rating_curve,\n terminal=terminal,\n discrete_control=discrete_control,\n starttime=\"2020-01-01 00:00:00\",\n endtime=\"2021-01-01 00:00:00\",\n)\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n<Axes: >\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_setpoint_with_minmax/ribasim.toml\")\n\nPosixPath('data/level_setpoint_with_minmax/ribasim.toml')\n\n\nNow run the model with level_setpoint_with_minmax/ribasim.toml. After running the model, read back the results:\n\nfrom matplotlib.dates import date2num\n\ndf_basin = pd.read_feather(datadir / \"level_setpoint_with_minmax/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\ndf_control = pd.read_feather(\n datadir / \"level_setpoint_with_minmax/results/control.arrow\"\n)\n\ny_min, y_max = ax.get_ybound()\nax.fill_between(df_control.time[:2], 2 * [y_min], 2 * [y_max], alpha=0.2, color=\"C0\")\nax.fill_between(df_control.time[2:4], 2 * [y_min], 2 * [y_max], alpha=0.2, color=\"C0\")\n\nax.set_xticks(\n date2num(df_control.time).tolist(),\n df_control.control_state.tolist(),\n rotation=50,\n)\n\nax.set_yticks(greater_than, [\"min\", \"setpoint\", \"max\"])\nax.set_ylabel(\"level\")\nplt.show()\n\n\n\n\nThe highlighted regions show where a pump is active.\nLet’s print an overview of what happened with control:\n\nmodel.print_discrete_control_record(\n datadir / \"level_setpoint_with_minmax/results/control.arrow\"\n)\n\n0. At 2020-01-01 00:00:00 the control node with ID 7 reached truth state TTT:\n For node ID 1 (Basin): level > 5.0\n For node ID 1 (Basin): level > 10.0\n For node ID 1 (Basin): level > 15.0\n\n This yielded control state \"out\":\n For node ID 2 (Pump): active = <NA>, flow_rate = 0.0, min_flow_rate = nan, max_flow_rate = nan\n For node ID 3 (Pump): active = <NA>, flow_rate = 0.002, min_flow_rate = nan, max_flow_rate = nan\n\n1. At 2020-02-08 19:02:21.861000 the control node with ID 7 reached truth state TFF:\n For node ID 1 (Basin): level > 5.0\n For node ID 1 (Basin): level < 10.0\n For node ID 1 (Basin): level < 15.0\n\n This yielded control state \"none\":\n For node ID 2 (Pump): active = <NA>, flow_rate = 0.0, min_flow_rate = nan, max_flow_rate = nan\n For node ID 3 (Pump): active = <NA>, flow_rate = 0.0, min_flow_rate = nan, max_flow_rate = nan\n\n2. At 2020-07-05 08:56:10.319000 the control node with ID 7 reached truth state FFF:\n For node ID 1 (Basin): level < 5.0\n For node ID 1 (Basin): level < 10.0\n For node ID 1 (Basin): level < 15.0\n\n This yielded control state \"in\":\n For node ID 2 (Pump): active = <NA>, flow_rate = 0.002, min_flow_rate = nan, max_flow_rate = nan\n For node ID 3 (Pump): active = <NA>, flow_rate = 0.0, min_flow_rate = nan, max_flow_rate = nan\n\n3. At 2020-08-11 06:05:15.592000 the control node with ID 7 reached truth state TTF:\n For node ID 1 (Basin): level > 5.0\n For node ID 1 (Basin): level > 10.0\n For node ID 1 (Basin): level < 15.0\n\n This yielded control state \"none\":\n For node ID 2 (Pump): active = <NA>, flow_rate = 0.0, min_flow_rate = nan, max_flow_rate = nan\n For node ID 3 (Pump): active = <NA>, flow_rate = 0.0, min_flow_rate = nan, max_flow_rate = nan\n\n\n\nNote that crossing direction specific truth states (containing “U”, “D”) are not present in this overview even though they are part of the control logic. This is because in the control logic for this model these truth states are only used to sustain control states, while the overview only shows changes in control states.\n\n\n4 Model with PID control\nSet up the nodes:\n\nxy = np.array(\n [\n (0.0, 0.0), # 1: FlowBoundary\n (1.0, 0.0), # 2: Basin\n (2.0, 0.5), # 3: Pump\n (3.0, 0.0), # 4: LevelBoundary\n (1.5, 1.0), # 5: PidControl\n (2.0, -0.5), # 6: outlet\n (1.5, -1.0), # 7: PidControl\n ]\n)\n\nnode_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n\nnode_type = [\n \"FlowBoundary\",\n \"Basin\",\n \"Pump\",\n \"LevelBoundary\",\n \"PidControl\",\n \"Outlet\",\n \"PidControl\",\n]\n\n# Make sure the feature id starts at 1: explicitly give an index.\nnode = ribasim.Node(\n df=gpd.GeoDataFrame(\n data={\"type\": node_type},\n index=pd.Index(np.arange(len(xy)) + 1, name=\"fid\"),\n geometry=node_xy,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the edges:\n\nfrom_id = np.array([1, 2, 3, 4, 6, 5, 7], dtype=np.int64)\nto_id = np.array([2, 3, 4, 6, 2, 3, 6], dtype=np.int64)\n\nlines = node.geometry_from_connectivity(from_id, to_id)\nedge = ribasim.Edge(\n df=gpd.GeoDataFrame(\n data={\n \"from_node_id\": from_id,\n \"to_node_id\": to_id,\n \"edge_type\": 5 * [\"flow\"] + 2 * [\"control\"],\n },\n geometry=lines,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the basins:\n\nprofile = pd.DataFrame(\n data={\"node_id\": [2, 2], \"level\": [0.0, 1.0], \"area\": [1000.0, 1000.0]}\n)\n\nstatic = pd.DataFrame(\n data={\n \"node_id\": [2],\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\nstate = pd.DataFrame(\n data={\n \"node_id\": [2],\n \"level\": [6.0],\n }\n)\n\nbasin = ribasim.Basin(profile=profile, static=static, state=state)\n\nSetup the pump:\n\npump = ribasim.Pump(\n static=pd.DataFrame(\n data={\n \"node_id\": [3],\n \"flow_rate\": [0.0], # Will be overwritten by PID controller\n }\n )\n)\n\nSetup the outlet:\n\noutlet = ribasim.Outlet(\n static=pd.DataFrame(\n data={\n \"node_id\": [6],\n \"flow_rate\": [0.0], # Will be overwritten by PID controller\n }\n )\n)\n\nSetup flow boundary:\n\nflow_boundary = ribasim.FlowBoundary(\n static=pd.DataFrame(data={\"node_id\": [1], \"flow_rate\": [1e-3]})\n)\n\nSetup flow boundary:\n\nlevel_boundary = ribasim.LevelBoundary(\n static=pd.DataFrame(\n data={\n \"node_id\": [4],\n \"level\": [1.0], # Not relevant\n }\n )\n)\n\nSetup PID control:\n\npid_control = ribasim.PidControl(\n time=pd.DataFrame(\n data={\n \"node_id\": 4 * [5, 7],\n \"time\": [\n \"2020-01-01 00:00:00\",\n \"2020-01-01 00:00:00\",\n \"2020-05-01 00:00:00\",\n \"2020-05-01 00:00:00\",\n \"2020-07-01 00:00:00\",\n \"2020-07-01 00:00:00\",\n \"2020-12-01 00:00:00\",\n \"2020-12-01 00:00:00\",\n ],\n \"listen_node_id\": 4 * [2, 2],\n \"target\": [5.0, 5.0, 5.0, 5.0, 7.5, 7.5, 7.5, 7.5],\n \"proportional\": 4 * [-1e-3, 1e-3],\n \"integral\": 4 * [-1e-7, 1e-7],\n \"derivative\": 4 * [0.0, 0.0],\n }\n )\n)\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 a model:\n\nmodel = ribasim.Model(\n network=ribasim.Network(\n node=node,\n edge=edge,\n ),\n basin=basin,\n flow_boundary=flow_boundary,\n level_boundary=level_boundary,\n pump=pump,\n outlet=outlet,\n pid_control=pid_control,\n starttime=\"2020-01-01 00:00:00\",\n endtime=\"2020-12-01 00:00:00\",\n)\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n<Axes: >\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\ntarget_levels = model.pid_control.time.df.target.to_numpy()[:4]\ntimes = date2num(model.pid_control.time.df.time)[:4]\nax.plot(times, target_levels, color=\"k\", ls=\":\", label=\"target level\")\npass\n\n\n\n\n\n\n5 Model with allocation\nSetup the nodes:\n\nxy = np.array(\n [\n (0.0, 0.0), # 1: FlowBoundary\n (1.0, 0.0), # 2: Basin\n (1.0, 1.0), # 3: User\n (2.0, 0.0), # 4: LinearResistance\n (3.0, 0.0), # 5: Basin\n (3.0, 1.0), # 6: User\n (4.0, 0.0), # 7: TabulatedRatingCurve\n (4.5, 0.0), # 8: FractionalFlow\n (4.5, 0.5), # 9: FractionalFlow\n (5.0, 0.0), # 10: Terminal\n (4.5, 0.25), # 11: DiscreteControl\n (4.5, 1.0), # 12: Basin\n (5.0, 1.0), # 13: User\n ]\n)\nnode_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])\n\nnode_type = [\n \"FlowBoundary\",\n \"Basin\",\n \"User\",\n \"LinearResistance\",\n \"Basin\",\n \"User\",\n \"TabulatedRatingCurve\",\n \"FractionalFlow\",\n \"FractionalFlow\",\n \"Terminal\",\n \"DiscreteControl\",\n \"Basin\",\n \"User\",\n]\n\n# All nodes belong to allocation network id 1\nnode = ribasim.Node(\n df=gpd.GeoDataFrame(\n data={\"type\": node_type, \"allocation_network_id\": 1},\n index=pd.Index(np.arange(len(xy)) + 1, name=\"fid\"),\n geometry=node_xy,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the edges:\n\nfrom_id = np.array(\n [1, 2, 2, 4, 5, 5, 7, 3, 6, 7, 8, 9, 12, 13, 11, 11],\n dtype=np.int64,\n)\nto_id = np.array(\n [2, 3, 4, 5, 6, 7, 8, 2, 5, 9, 10, 12, 13, 10, 8, 9],\n dtype=np.int64,\n)\n# Denote the first edge, 1 => 2, as a source edge for\n# allocation network 1\nallocation_network_id = len(from_id) * [None]\nallocation_network_id[0] = 1\nlines = node.geometry_from_connectivity(from_id, to_id)\nedge = ribasim.Edge(\n df=gpd.GeoDataFrame(\n data={\n \"from_node_id\": from_id,\n \"to_node_id\": to_id,\n \"edge_type\": (len(from_id) - 2) * [\"flow\"] + 2 * [\"control\"],\n \"allocation_network_id\": allocation_network_id,\n },\n geometry=lines,\n crs=\"EPSG:28992\",\n )\n)\n\nSetup the basins:\n\nprofile = pd.DataFrame(\n data={\n \"node_id\": [2, 2, 5, 5, 12, 12],\n \"area\": 300_000.0,\n \"level\": 3 * [0.0, 1.0],\n }\n)\n\nstatic = pd.DataFrame(\n data={\n \"node_id\": [2, 5, 12],\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\nstate = pd.DataFrame(data={\"node_id\": [2, 5, 12], \"level\": 1.0})\n\nbasin = ribasim.Basin(profile=profile, static=static, state=state)\n\nSetup the flow boundary:\n\nflow_boundary = ribasim.FlowBoundary(\n static=pd.DataFrame(\n data={\n \"node_id\": [1],\n \"flow_rate\": 2.0,\n }\n )\n)\n\nSetup the linear resistance:\n\nlinear_resistance = ribasim.LinearResistance(\n static=pd.DataFrame(\n data={\n \"node_id\": [4],\n \"resistance\": 0.06,\n }\n )\n)\n\nSetup the tabulated rating curve:\n\ntabulated_rating_curve = ribasim.TabulatedRatingCurve(\n static=pd.DataFrame(\n data={\n \"node_id\": 7,\n \"level\": [0.0, 0.5, 1.0],\n \"flow_rate\": [0.0, 0.0, 2.0],\n }\n )\n)\n\nSetup the fractional flow:\n\nfractional_flow = ribasim.FractionalFlow(\n static=pd.DataFrame(\n data={\n \"node_id\": [8, 8, 9, 9],\n \"fraction\": [0.6, 0.9, 0.4, 0.1],\n \"control_state\": [\"divert\", \"close\", \"divert\", \"close\"],\n }\n )\n)\n\nSetup the terminal:\n\nterminal = ribasim.Terminal(\n static=pd.DataFrame(\n data={\n \"node_id\": [10],\n }\n )\n)\n\nSetup the discrete control:\n\ncondition = pd.DataFrame(\n data={\n \"node_id\": [11],\n \"listen_feature_id\": 5,\n \"variable\": \"level\",\n \"greater_than\": 0.52,\n }\n)\n\nlogic = pd.DataFrame(\n data={\n \"node_id\": 11,\n \"truth_state\": [\"T\", \"F\"],\n \"control_state\": [\"divert\", \"close\"],\n }\n)\n\ndiscrete_control = ribasim.DiscreteControl(condition=condition, logic=logic)\n\nSetup the users:\n\nuser = ribasim.User(\n static=pd.DataFrame(\n data={\n \"node_id\": [6, 13],\n \"demand\": [1.5, 1.0],\n \"return_factor\": 0.0,\n \"min_level\": -1.0,\n \"priority\": [1, 3],\n }\n ),\n time=pd.DataFrame(\n data={\n \"node_id\": [3, 3, 3, 3],\n \"demand\": [0.0, 1.0, 1.2, 1.2],\n \"priority\": [1, 1, 2, 2],\n \"return_factor\": 0.0,\n \"min_level\": -1.0,\n \"time\": 2 * [\"2020-01-01 00:00:00\", \"2020-01-20 00:00:00\"],\n }\n ),\n)\n\nSetup the allocation:\n\nallocation = ribasim.Allocation(use_allocation=True, timestep=86400)\n\nSetup a model:\n\nmodel = ribasim.Model(\n network=ribasim.Network(\n node=node,\n edge=edge,\n ),\n basin=basin,\n flow_boundary=flow_boundary,\n linear_resistance=linear_resistance,\n tabulated_rating_curve=tabulated_rating_curve,\n terminal=terminal,\n user=user,\n discrete_control=discrete_control,\n fractional_flow=fractional_flow,\n allocation=allocation,\n starttime=\"2020-01-01 00:00:00\",\n endtime=\"2020-01-20 00:00:00\",\n)\n\nLet’s take a look at the model:\n\nmodel.plot()\n\n<Axes: >\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=[\"user_node_id\", \"priority\"],\n values=[\"demand\", \"allocated\", \"abstracted\"],\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))\nfig.suptitle(f\"Objective type: {model.allocation.objective_type}\", fontweight=\"bold\")\n\ndf_allocation_wide[\"demand\"].plot(ax=axs[0], ls=\":\")\ndf_allocation_wide[\"allocated\"].plot(ax=axs[1], ls=\"--\")\ndf_allocation_wide[\"abstracted\"].plot(ax=axs[2])\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\nSome things to note about this plot:\n\nAbstraction behaves somewhat erratically at the start of the simulation. This is because allocation is based on flows computed in the physical layer, and at the start of the simulation these are not known yet.\nAlthough there is a plotted line for abstraction per priority, abstraction is actually 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]')"
},
{
"objectID": "src/index.html",
@@ -606,7 +606,7 @@
"href": "core/allocation.html#example",
"title": "Allocation",
"section": "5.1 Example",
- "text": "5.1 Example\nThe following is an example of an optimization problem for the example shown here:\n\n\nCode\nusing Ribasim\nusing SQLite\n\ntoml_path = normpath(@__DIR__, \"../../generated_testmodels/allocation_example/ribasim.toml\")\np = Ribasim.Model(toml_path).integrator.p\n\nallocation_model = p.allocation.allocation_models[1]\nt = 0.0\npriority_idx = 1\n\nRibasim.set_flow!(p.graph, Ribasim.NodeID(1), Ribasim.NodeID(2), 1.0)\n\nRibasim.adjust_source_capacities!(allocation_model, p, priority_idx)\nRibasim.adjust_edge_capacities!(allocation_model, p, priority_idx)\nRibasim.set_objective_priority!(allocation_model, p, t, priority_idx)\n\nprintln(p.allocation.allocation_models[1].problem)\n\n\nMin 0.05 F[(#1, #2)] + 0.05 F[(#2, #3)] + 0.05 F[(#12, #13)] + 0.05 F[(#2, #5)] + 0.05 F[(#5, #7)] + 0.05 F[(#7, #12)] + 0.05 F[(#13, #10)] + 0.05 F[(#5, #6)] + 0.05 F[(#7, #10)] + 0.05 F[(#5, #2)] + F_abs[#13] + F_abs[#6] + F_abs[#3]\nSubject to\n source[(#1, #2)] : F[(#1, #2)] ≤ 1\n flow_conservation[#12] : F[(#12, #13)] - F[(#7, #12)] ≤ 0\n flow_conservation[#2] : -F[(#1, #2)] + F[(#2, #3)] + F[(#2, #5)] - F[(#5, #2)] ≤ 0\n flow_conservation[#5] : -F[(#2, #5)] + F[(#5, #7)] + F[(#5, #6)] + F[(#5, #2)] ≤ 0\n return_flow[#13] : F[(#13, #10)] ≤ 0\n fractional_flow[(#7, #12)] : -0.4 F[(#5, #7)] + F[(#7, #12)] ≤ 0\n F[(#1, #2)] ≥ 0\n F[(#2, #3)] ≥ 0\n F[(#12, #13)] ≥ 0\n F[(#2, #5)] ≥ 0\n F[(#5, #7)] ≥ 0\n F[(#7, #12)] ≥ 0\n F[(#13, #10)] ≥ 0\n F[(#5, #6)] ≥ 0\n F[(#7, #10)] ≥ 0\n F[(#5, #2)] ≥ 0\n abs_positive[#13] : -F[(#12, #13)] + F_abs[#13] ≥ 0\n abs_positive[#6] : -F[(#5, #6)] + F_abs[#6] ≥ -1.5\n abs_positive[#3] : -F[(#2, #3)] + F_abs[#3] ≥ 0\n abs_negative[#13] : F[(#12, #13)] + F_abs[#13] ≥ 0\n abs_negative[#6] : F[(#5, #6)] + F_abs[#6] ≥ 1.5\n abs_negative[#3] : F[(#2, #3)] + F_abs[#3] ≥ 0"
+ "text": "5.1 Example\nThe following is an example of an optimization problem for the example shown here:\n\n\nCode\nusing Ribasim\nusing SQLite\n\ntoml_path = normpath(@__DIR__, \"../../generated_testmodels/allocation_example/ribasim.toml\")\np = Ribasim.Model(toml_path).integrator.p\n\nallocation_model = p.allocation.allocation_models[1]\nt = 0.0\npriority_idx = 1\n\nRibasim.set_flow!(p.graph, Ribasim.NodeID(1), Ribasim.NodeID(2), 1.0)\n\nRibasim.adjust_source_capacities!(allocation_model, p, priority_idx)\nRibasim.adjust_edge_capacities!(allocation_model, p, priority_idx)\nRibasim.set_objective_priority!(allocation_model, p, t, priority_idx)\n\nprintln(p.allocation.allocation_models[1].problem)\n\n\nMin 0.05 F[(#2, #3)] + 0.05 F[(#7, #10)] + 0.05 F[(#5, #7)] + 0.05 F[(#1, #2)] + 0.05 F[(#2, #5)] + 0.05 F[(#7, #12)] + 0.05 F[(#5, #2)] + 0.05 F[(#13, #10)] + 0.05 F[(#5, #6)] + 0.05 F[(#12, #13)] + F_abs[#6] + F_abs[#13] + F_abs[#3]\nSubject to\n source[(#1, #2)] : F[(#1, #2)] ≤ 1\n flow_conservation[#5] : F[(#5, #7)] - F[(#2, #5)] + F[(#5, #2)] + F[(#5, #6)] ≤ 0\n flow_conservation[#12] : -F[(#7, #12)] + F[(#12, #13)] ≤ 0\n flow_conservation[#2] : F[(#2, #3)] - F[(#1, #2)] + F[(#2, #5)] - F[(#5, #2)] ≤ 0\n return_flow[#13] : F[(#13, #10)] ≤ 0\n fractional_flow[(#7, #12)] : -0.4 F[(#5, #7)] + F[(#7, #12)] ≤ 0\n F[(#2, #3)] ≥ 0\n F[(#7, #10)] ≥ 0\n F[(#5, #7)] ≥ 0\n F[(#1, #2)] ≥ 0\n F[(#2, #5)] ≥ 0\n F[(#7, #12)] ≥ 0\n F[(#5, #2)] ≥ 0\n F[(#13, #10)] ≥ 0\n F[(#5, #6)] ≥ 0\n F[(#12, #13)] ≥ 0\n abs_positive[#6] : -F[(#5, #6)] + F_abs[#6] ≥ -1.5\n abs_positive[#13] : -F[(#12, #13)] + F_abs[#13] ≥ 0\n abs_positive[#3] : -F[(#2, #3)] + F_abs[#3] ≥ 0\n abs_negative[#6] : F[(#5, #6)] + F_abs[#6] ≥ 1.5\n abs_negative[#13] : F[(#12, #13)] + F_abs[#13] ≥ 0\n abs_negative[#3] : F[(#2, #3)] + F_abs[#3] ≥ 0"
},
{
"objectID": "core/usage.html",
@@ -830,7 +830,7 @@
"href": "core/equations.html#precipitation",
"title": "Equations",
"section": "2.2 Precipitation",
- "text": "2.2 Precipitation\nThe precipitation term is given by\n\\[\n Q_P = P \\cdot A(u).\n\\tag{2}\\]\nHere \\(P = P(t)\\) is the precipitation rate and \\(A\\) is the wetted area. \\(A\\) is a function of the storage \\(u = u(t)\\): as the volume of water changes, the area of the free water surface may change as well, depending on the slopes of the surface waters."
+ "text": "2.2 Precipitation\nThe precipitation term is given by\n\\[\n Q_P = P \\cdot A.\n\\tag{2}\\]\nHere \\(P = P(t)\\) is the precipitation rate and \\(A\\) is the maximum area given in the Basin / profile table. Precipitation in the Basin area is assumed to be directly added to the Basin storage. The modeler needs to ensure all precipitation enters the model, and there is no overlap in the maximum profile areas, else extra water is created. If a part of the catchment is not in any Basin profile, the modeler has to verify that water source is not forgotten. It can for instance be converted to a flow rate and added to a Basin as a FlowBoundary."
},
{
"objectID": "core/equations.html#evaporation",
@@ -844,28 +844,28 @@
"href": "core/equations.html#infiltration-and-drainage",
"title": "Equations",
"section": "2.4 Infiltration and Drainage",
- "text": "2.4 Infiltration and Drainage\nInfiltration is provided as a lump sum for the basin. If Ribasim is coupled with MODFLOW 6, the infiltration is computed as the sum of all positive flows of the MODFLOW 6 boundary conditions in the basin:\n\\[\n Q_\\text{inf} = \\sum_{i=1}^{n} \\sum_{j=1}^{m} \\max(Q_{\\mathrm{mf6}_{i,j}}, 0.0)\n\\] {#eq-inf}.\nWhere \\(i\\) is the index of the boundary condition, \\(j\\) the MODFLOW 6 cell index, \\(n\\) the number of boundary conditions, and \\(m\\) the number of MODFLOW 6 cells in the basin. \\(Q_{\\mathrm{mf6}_{i,j}}\\) is the flow computed by MODFLOW 6 for cell \\(j\\) for boundary condition \\(i\\).\nDrainage is a lump sump for the basin, and consists of the sum of the absolute value of all negative flows of the MODFLOW 6 boundary conditions in the basin.\n\\[\n Q_\\text{drn} = \\sum_{i=1}^{n} \\sum_{j=1}^{m} \\left| \\min(Q_{\\mathrm{mf6}_{i,j}}, 0.0) \\right|\n\\tag{4}\\]\nThe interaction with MODFLOW 6 boundary conditions is explained in greater detail in the the MODFLOW coupling section of the documentation."
+ "text": "2.4 Infiltration and Drainage\nInfiltration is provided as a lump sum for the basin. If Ribasim is coupled with MODFLOW 6, the infiltration is computed as the sum of all positive flows of the MODFLOW 6 boundary conditions in the basin:\n\\[\n Q_\\text{inf} = \\sum_{i=1}^{n} \\sum_{j=1}^{m} \\max(Q_{\\mathrm{mf6}_{i,j}}, 0.0)\n\\tag{4}\\]\nWhere \\(i\\) is the index of the boundary condition, \\(j\\) the MODFLOW 6 cell index, \\(n\\) the number of boundary conditions, and \\(m\\) the number of MODFLOW 6 cells in the basin. \\(Q_{\\mathrm{mf6}_{i,j}}\\) is the flow computed by MODFLOW 6 for cell \\(j\\) for boundary condition \\(i\\).\nDrainage is a lump sump for the basin, and consists of the sum of the absolute value of all negative flows of the MODFLOW 6 boundary conditions in the basin.\n\\[\n Q_\\text{drn} = \\sum_{i=1}^{n} \\sum_{j=1}^{m} \\left| \\min(Q_{\\mathrm{mf6}_{i,j}}, 0.0) \\right|\n\\tag{5}\\]\nThe interaction with MODFLOW 6 boundary conditions is explained in greater detail in the the MODFLOW coupling section of the documentation."
},
{
"objectID": "core/equations.html#upstream-and-downstream-flow",
"href": "core/equations.html#upstream-and-downstream-flow",
"title": "Equations",
"section": "2.5 Upstream and downstream flow",
- "text": "2.5 Upstream and downstream flow\nRibasim’s basins can be connected to each other, and each basin expects an explicit connection. These connections are currently available for inter-basin flows:\n\n\nPump\nTabulatedRatingCurve\nLinearResistance\nManningResistance\n\nThe flow direction of the basin is not pre-determined: flow directions may freely reverse, provided the connection allows it. Currently, a LinearResistance allows bidirectional flow, but the\nAdditionally, three additional “connections” area available for the “outmost” basins (external nodes) in a network.\n\nTerminal\nLevelBoundary\nFlowBoundary\n\n\n2.5.1 Pump\nThe behaviour of pumps is very straight forward if these nodes are not PID controlled. Their flow is given by a fixed flow rate \\(q\\), multiplied by a reduction factor: \\[\nQ_\\text{pump} = \\phi(u; 10.0)q\n\\]\nHere \\(u\\) is the storage of the upstream basin. The reduction factor \\(\\phi\\) makes sure that the flow of the pump goes smootly to \\(0\\) as the upstream basin dries out.\n\n\n2.5.2 Outlet\nThe outlet is very similar to the pump, but it has a few extra reduction factors for physical constraints: \\[\nQ_\\text{outlet} = \\phi(u_a; 10.0)\\phi(\\Delta h; 0.1) \\phi(h_a-h_\\text{min};0.1)q.\n\\] The subscript \\(a\\) denotes the upstream node and \\(b\\) the downstream node. The first reduction factor is equivalent to the one for the pump. The second one makes sure that the outlet flow goes to zero as the head difference \\(\\Delta h = h_a - h_b\\) goes to zero. The last one makes sure that the outlet only produces flow when the upstream level is above the minimum chrest level \\(h_\\text{min}\\).\nNot all node types upstream or downstream of the outlet have a defined level. If this is the case, and therefore the reduction factor cannot be computed, it is defined to be \\(1.0\\).\n\n\n2.5.3 TabulatedRatingCurve\nThe Tabulated Rating Curve is a tabulation of a basin’s discharge behavior. It describes a piecewise linear relationship between the basin’s level and its discharge. It can be understood as an empirical description of a basin’s properties. This can include an outlet, but also the lumped hydraulic behavior of the upstream channels.\n\n\n\n\n\n\nNote\n\n\n\nCurrently, the discharge relies only on the basin’s level; it could also use the volume of both connected basins to simulate backwater effects, submersion of outlets, or even reversal of flows for high precipitation events.\n\n\n\n\n2.5.4 LinearResistance\nA LinearResistance connects two basins together. The flow between the two basins is determined by a linear relationship:\n\\[\n Q = \\frac{h_a - h_b}{R}\n\\tag{5}\\]\nHere \\(h_a\\) is the water level in the first basin, \\(h_b\\) is the water level in the second basin, and \\(R\\) is the resistance of the link. A LinearResistance makes no assumptions about the direction of the flow: water flows from high to low.\n\n\n2.5.5 Terminal\nThis only allows outflow from a basin into a terminal node.\n\n\n2.5.6 LevelBoundary\nThis can be connected to a basin via a LinearResistance. This boundary node will then exchange water with the basin based on the difference in water level between the two.\n\n\n2.5.7 FlowBoundary\nThis can be connected directly to a basin and prescribes the flow to or from that basin. We require that the edge connecting the flow boundary to the basin should point towards the basin, so that positive flow corresponds to water being added to the model.\n\n\n2.5.8 Manning connection\nRibasim is capable of simulating steady flow between basins through a reach described by a trapezoidal profile and a Manning roughness coefficient.\nWe describe the discharge from basin \\(a\\) to basin \\(b\\) solely as a function of the water levels in \\(a\\) and \\(b\\).\n\\[\nQ = f(h_a, h_b)\n\\]\nwhere:\n\nThe subscripts \\(a,b\\) denote basins\n\\(h\\) is the hydraulic head, or water level\n\nThe energy equation for open channel flow is:\n\\[\nH = h + \\frac{v^2}{2g}\n\\]\nWhere\n\n\\(H\\) is total head\n\\(v\\) is average water velocity\n\\(g\\) is gravitational acceleration\n\nThe discharge \\(Q\\) is defined as:\n\\[\nQ = Av\n\\]\nwhere \\(A\\) is cross-sectional area.\nWe use conservation of energy to relate the total head at \\(a\\) to \\(b\\), with \\(H_a > H_b\\) as follows:\n\\[\nH_a = H_b + h_{\\text{loss}}\n\\]\nOr:\n\\[\nh_a + \\frac{v_a^2}{2g} = h_b + \\frac{v_b^2}{2g} + h_{\\text{loss}}\n\\]\nWhere \\(v\\) is the average water velocity. \\(h_{\\text{loss}}\\) is a combination of friction and contraction/expansion losses:\n\\[\nh_{\\text{loss}} = S_f L + \\frac{C}{2g} \\left(v_b^2 - v_a^2\\right)\n\\]\nWhere:\n\n\\(L\\) is the reach length\n\\(S_f\\) is the representative friction slope\n\\(C\\) is the expansion or contraction coefficient, \\(0 \\le C \\le1\\)\n\nWe assume velocity differences in a connection are negligible (\\(v_a = v_b\\)):\n\\[\nh_a = h_b + S_f L\n\\]\nFriction losses are computed with the Gauckler-Manning formula:\n\\[\nQ = \\frac{A}{n} R_h^\\frac{2}{3} \\sqrt{S_f}\n\\]\nWhere:\n\n\\(A\\) is the representative area.\n\\(R_h\\) is the representative wetted radius.\n\\(S_f\\) is the representative friction slope.\n\\(n\\) is Manning’s roughness coefficient.\n\nWe can rewrite to express \\(S_f\\) in terms of Q:\n\\[\nS_f = Q^2 \\frac{n^2}{A^2 R_h^{4/3}}\n\\]\nNo water is added or removed in a connection:\n\\[\nQ_a = Q_b = Q\n\\]\nSubstituting:\n\\[\nh_a = h_b + Q^2 \\frac{n^2}{A^2 R_h^{4/3}} L\n\\]\nWe can then express \\(Q\\) as a function of head difference \\(\\Delta h\\):\n\\[\nQ = \\textrm{sign}(\\Delta h) \\frac{A}{n} R_h^{2/3}\\sqrt{\\frac{|\\Delta h|}{L} }\n\\]\nThe \\(\\textrm{sign}(\\Delta h)\\) term causes the direction of the flow to reverse if the head in basin \\(b\\) is larger than in basin \\(a\\).\nThis expression however leads to problems in simulation since the derivative of \\(Q\\) with respect to \\(\\Delta h\\) tends to \\(\\pm \\infty\\) as \\(\\Delta h\\) tends to 0. Therefore we use the slightly modified expression\n\\[\nQ = \\textrm{sign}(\\Delta h) \\frac{A}{n} R_h^{2/3}\\sqrt{\\frac{\\Delta h}{L} s(\\Delta h)}\n\\]\nto smooth out this problem. Here \\(s(x) = \\frac{2}{\\pi}\\arctan{1000x}\\) can be thought of as a smooth approximation of the sign function.\n\n\n\n\n\n\nNote\n\n\n\nThe computation of \\(S_f\\) is not exact: we base it on a representative area and hydraulic radius, rather than integrating \\(S_f\\) along the length of a reach. Direct analytic solutions exist for e.g. parabolic profiles (Tolkmitt), but other profiles requires relatively complicated approaches (such as approximating the profile with a polynomial).\nWe use the average value of the cross-sectional area, the average value of the water depth, and the average value of the hydraulic radius to compute a friction slope. The size of the resulting error will depend on the water depth difference between the upstream and downstream basin.\n\n\nThe cross sectional area for a trapezoidal or rectangular profile:\n\\[\nA = w d + \\frac{\\Delta y}{\\Delta z} d^2\n\\]\nWhere\n\n\\(w\\) is the width at \\(d = 0\\) (A triangular profile has \\(w = 0\\))\n\\(\\frac{\\Delta y}{\\Delta z}\\) is the slope of the profile expressed as the horizontal length for one unit in the vertical (A slope of 45 degrees has \\(\\frac{\\Delta y}{\\Delta z} = 1\\); a rectangular profile 0).\n\nAccordingly, the wetted perimeter is:\n\\[\nB = w + 2 d \\sqrt{\\left(\\frac{\\Delta y}{\\Delta z}\\right)^2 + 1}\n\\]"
+ "text": "2.5 Upstream and downstream flow\nRibasim’s basins can be connected to each other, and each basin expects an explicit connection. These connections are currently available for inter-basin flows:\n\n\nPump\nTabulatedRatingCurve\nLinearResistance\nManningResistance\n\nThe flow direction of the basin is not pre-determined: flow directions may freely reverse, provided the connection allows it. Currently, a LinearResistance allows bidirectional flow, but the\nAdditionally, three additional “connections” area available for the “outmost” basins (external nodes) in a network.\n\nTerminal\nLevelBoundary\nFlowBoundary\n\n\n2.5.1 Pump\nThe behaviour of pumps is very straight forward if these nodes are not PID controlled. Their flow is given by a fixed flow rate \\(q\\), multiplied by a reduction factor: \\[\nQ_\\text{pump} = \\phi(u; 10.0)q\n\\]\nHere \\(u\\) is the storage of the upstream basin. The reduction factor \\(\\phi\\) makes sure that the flow of the pump goes smootly to \\(0\\) as the upstream basin dries out.\n\n\n2.5.2 Outlet\nThe outlet is very similar to the pump, but it has a few extra reduction factors for physical constraints: \\[\nQ_\\text{outlet} = \\phi(u_a; 10.0)\\phi(\\Delta h; 0.1) \\phi(h_a-h_\\text{min};0.1)q.\n\\] The subscript \\(a\\) denotes the upstream node and \\(b\\) the downstream node. The first reduction factor is equivalent to the one for the pump. The second one makes sure that the outlet flow goes to zero as the head difference \\(\\Delta h = h_a - h_b\\) goes to zero. The last one makes sure that the outlet only produces flow when the upstream level is above the minimum chrest level \\(h_\\text{min}\\).\nNot all node types upstream or downstream of the outlet have a defined level. If this is the case, and therefore the reduction factor cannot be computed, it is defined to be \\(1.0\\).\n\n\n2.5.3 TabulatedRatingCurve\nThe Tabulated Rating Curve is a tabulation of a basin’s discharge behavior. It describes a piecewise linear relationship between the basin’s level and its discharge. It can be understood as an empirical description of a basin’s properties. This can include an outlet, but also the lumped hydraulic behavior of the upstream channels.\n\n\n\n\n\n\nNote\n\n\n\nCurrently, the discharge relies only on the basin’s level; it could also use the volume of both connected basins to simulate backwater effects, submersion of outlets, or even reversal of flows for high precipitation events.\n\n\n\n\n2.5.4 LinearResistance\nA LinearResistance connects two basins together. The flow between the two basins is determined by a linear relationship:\n\\[\n Q = \\frac{h_a - h_b}{R}\n\\tag{6}\\]\nHere \\(h_a\\) is the water level in the first basin, \\(h_b\\) is the water level in the second basin, and \\(R\\) is the resistance of the link. A LinearResistance makes no assumptions about the direction of the flow: water flows from high to low.\n\n\n2.5.5 Terminal\nThis only allows outflow from a basin into a terminal node.\n\n\n2.5.6 LevelBoundary\nThis can be connected to a basin via a LinearResistance. This boundary node will then exchange water with the basin based on the difference in water level between the two.\n\n\n2.5.7 FlowBoundary\nThis can be connected directly to a basin and prescribes the flow to or from that basin. We require that the edge connecting the flow boundary to the basin should point towards the basin, so that positive flow corresponds to water being added to the model.\n\n\n2.5.8 Manning connection\nRibasim is capable of simulating steady flow between basins through a reach described by a trapezoidal profile and a Manning roughness coefficient.\nWe describe the discharge from basin \\(a\\) to basin \\(b\\) solely as a function of the water levels in \\(a\\) and \\(b\\).\n\\[\nQ = f(h_a, h_b)\n\\]\nwhere:\n\nThe subscripts \\(a,b\\) denote basins\n\\(h\\) is the hydraulic head, or water level\n\nThe energy equation for open channel flow is:\n\\[\nH = h + \\frac{v^2}{2g}\n\\]\nWhere\n\n\\(H\\) is total head\n\\(v\\) is average water velocity\n\\(g\\) is gravitational acceleration\n\nThe discharge \\(Q\\) is defined as:\n\\[\nQ = Av\n\\]\nwhere \\(A\\) is cross-sectional area.\nWe use conservation of energy to relate the total head at \\(a\\) to \\(b\\), with \\(H_a > H_b\\) as follows:\n\\[\nH_a = H_b + h_{\\text{loss}}\n\\]\nOr:\n\\[\nh_a + \\frac{v_a^2}{2g} = h_b + \\frac{v_b^2}{2g} + h_{\\text{loss}}\n\\]\nWhere \\(v\\) is the average water velocity. \\(h_{\\text{loss}}\\) is a combination of friction and contraction/expansion losses:\n\\[\nh_{\\text{loss}} = S_f L + \\frac{C}{2g} \\left(v_b^2 - v_a^2\\right)\n\\]\nWhere:\n\n\\(L\\) is the reach length\n\\(S_f\\) is the representative friction slope\n\\(C\\) is the expansion or contraction coefficient, \\(0 \\le C \\le1\\)\n\nWe assume velocity differences in a connection are negligible (\\(v_a = v_b\\)):\n\\[\nh_a = h_b + S_f L\n\\]\nFriction losses are computed with the Gauckler-Manning formula:\n\\[\nQ = \\frac{A}{n} R_h^\\frac{2}{3} \\sqrt{S_f}\n\\]\nWhere:\n\n\\(A\\) is the representative area.\n\\(R_h\\) is the representative wetted radius.\n\\(S_f\\) is the representative friction slope.\n\\(n\\) is Manning’s roughness coefficient.\n\nWe can rewrite to express \\(S_f\\) in terms of Q:\n\\[\nS_f = Q^2 \\frac{n^2}{A^2 R_h^{4/3}}\n\\]\nNo water is added or removed in a connection:\n\\[\nQ_a = Q_b = Q\n\\]\nSubstituting:\n\\[\nh_a = h_b + Q^2 \\frac{n^2}{A^2 R_h^{4/3}} L\n\\]\nWe can then express \\(Q\\) as a function of head difference \\(\\Delta h\\):\n\\[\nQ = \\textrm{sign}(\\Delta h) \\frac{A}{n} R_h^{2/3}\\sqrt{\\frac{|\\Delta h|}{L} }\n\\]\nThe \\(\\textrm{sign}(\\Delta h)\\) term causes the direction of the flow to reverse if the head in basin \\(b\\) is larger than in basin \\(a\\).\nThis expression however leads to problems in simulation since the derivative of \\(Q\\) with respect to \\(\\Delta h\\) tends to \\(\\pm \\infty\\) as \\(\\Delta h\\) tends to 0. Therefore we use the slightly modified expression\n\\[\nQ = \\textrm{sign}(\\Delta h) \\frac{A}{n} R_h^{2/3}\\sqrt{\\frac{\\Delta h}{L} s(\\Delta h)}\n\\]\nto smooth out this problem. Here \\(s(x) = \\frac{2}{\\pi}\\arctan{1000x}\\) can be thought of as a smooth approximation of the sign function.\n\n\n\n\n\n\nNote\n\n\n\nThe computation of \\(S_f\\) is not exact: we base it on a representative area and hydraulic radius, rather than integrating \\(S_f\\) along the length of a reach. Direct analytic solutions exist for e.g. parabolic profiles (Tolkmitt), but other profiles requires relatively complicated approaches (such as approximating the profile with a polynomial).\nWe use the average value of the cross-sectional area, the average value of the water depth, and the average value of the hydraulic radius to compute a friction slope. The size of the resulting error will depend on the water depth difference between the upstream and downstream basin.\n\n\nThe cross sectional area for a trapezoidal or rectangular profile:\n\\[\nA = w d + \\frac{\\Delta y}{\\Delta z} d^2\n\\]\nWhere\n\n\\(w\\) is the width at \\(d = 0\\) (A triangular profile has \\(w = 0\\))\n\\(\\frac{\\Delta y}{\\Delta z}\\) is the slope of the profile expressed as the horizontal length for one unit in the vertical (A slope of 45 degrees has \\(\\frac{\\Delta y}{\\Delta z} = 1\\); a rectangular profile 0).\n\nAccordingly, the wetted perimeter is:\n\\[\nB = w + 2 d \\sqrt{\\left(\\frac{\\Delta y}{\\Delta z}\\right)^2 + 1}\n\\]"
},
{
"objectID": "core/equations.html#the-derivative-term",
"href": "core/equations.html#the-derivative-term",
"title": "Equations",
"section": "4.1 The derivative term",
- "text": "4.1 The derivative term\nWhen \\(K_d \\ne 0\\) this adds a level of complexity. We can see this by looking at the error derivative more closely: \\[\n\\frac{\\text{d}e}{\\text{d}t} = \\frac{\\text{d}\\text{SP}}{\\text{d}t} - \\frac{1}{A(u_\\text{PID})}\\frac{\\text{d}u_\\text{PID}}{\\text{d}t},\n\\] where \\(A(u_\\text{PID})\\) is the area of the controlled basin as a function of the storage of the controlled basin \\(u_\\text{PID}\\). The complexity arises from the fact that \\(Q_\\text{PID}\\) is a contribution to \\(\\frac{\\text{d}u_\\text{PID}}{\\text{d}t} = f_\\text{PID}\\), which makes Equation 7 an implicit equation for \\(Q_\\text{PID}\\). We define\n\\[\nf_\\text{PID} = \\hat{f}_\\text{PID} \\pm Q_\\text{pump/outlet},\n\\]\nthat is, \\(\\hat{f}_\\text{PID}\\) is the right hand side of the ODE for the controlled basin storage state without the contribution of the PID controlled pump. The plus sign holds for an outlet and the minus sign for a pump, dictated by the way the pump and outlet connectivity to the controlled basin is enforced.\nUsing this, solving Equation 7 for \\(Q_\\text{PID}\\) yields \\[\nQ_\\text{pump/outlet} = \\text{clip}\\left(\\phi(u_\\text{us})\\frac{K_pe + K_iI + K_d \\left(\\frac{\\text{d}\\text{SP}}{\\text{d}t}-\\frac{\\hat{f}_\\text{PID}}{A(u_\\text{PID})}\\right)}{1\\pm\\phi(u_\\text{us})\\frac{K_d}{A(u_\\text{PID})}};Q_{\\min},Q_{\\max}\\right),\n\\] where the clipping is again done last. Note that to compute this, \\(\\hat{f}_\\text{PID}\\) has to be known first, meaning that the PID controlled pump/outlet flow rate has to be computed after all other contributions to the PID controlled basin’s storage are known."
+ "text": "4.1 The derivative term\nWhen \\(K_d \\ne 0\\) this adds a level of complexity. We can see this by looking at the error derivative more closely: \\[\n\\frac{\\text{d}e}{\\text{d}t} = \\frac{\\text{d}\\text{SP}}{\\text{d}t} - \\frac{1}{A(u_\\text{PID})}\\frac{\\text{d}u_\\text{PID}}{\\text{d}t},\n\\] where \\(A(u_\\text{PID})\\) is the area of the controlled basin as a function of the storage of the controlled basin \\(u_\\text{PID}\\). The complexity arises from the fact that \\(Q_\\text{PID}\\) is a contribution to \\(\\frac{\\text{d}u_\\text{PID}}{\\text{d}t} = f_\\text{PID}\\), which makes Equation 8 an implicit equation for \\(Q_\\text{PID}\\). We define\n\\[\nf_\\text{PID} = \\hat{f}_\\text{PID} \\pm Q_\\text{pump/outlet},\n\\]\nthat is, \\(\\hat{f}_\\text{PID}\\) is the right hand side of the ODE for the controlled basin storage state without the contribution of the PID controlled pump. The plus sign holds for an outlet and the minus sign for a pump, dictated by the way the pump and outlet connectivity to the controlled basin is enforced.\nUsing this, solving Equation 8 for \\(Q_\\text{PID}\\) yields \\[\nQ_\\text{pump/outlet} = \\text{clip}\\left(\\phi(u_\\text{us})\\frac{K_pe + K_iI + K_d \\left(\\frac{\\text{d}\\text{SP}}{\\text{d}t}-\\frac{\\hat{f}_\\text{PID}}{A(u_\\text{PID})}\\right)}{1\\pm\\phi(u_\\text{us})\\frac{K_d}{A(u_\\text{PID})}};Q_{\\min},Q_{\\max}\\right),\n\\] where the clipping is again done last. Note that to compute this, \\(\\hat{f}_\\text{PID}\\) has to be known first, meaning that the PID controlled pump/outlet flow rate has to be computed after all other contributions to the PID controlled basin’s storage are known."
},
{
"objectID": "core/equations.html#the-sign-of-the-parameters",
"href": "core/equations.html#the-sign-of-the-parameters",
"title": "Equations",
"section": "4.2 The sign of the parameters",
- "text": "4.2 The sign of the parameters\nNote by Equation 6 that the error is positive if the setpoint is larger than the basin level and negative if the setpoint is smaller than the basin level.\nWe enforce the convention that when a pump is controlled, its edge points away from the basin, and when an outlet is controlled, its edge points towards the basin, so that the main flow direction along these edges is positive. Therefore, positive flows of the pump and outlet have opposite effects on the basin, and thus the parameters \\(K_p,K_i,K_d\\) of the pump and outlet must have oppositive signs to achieve the same goal."
+ "text": "4.2 The sign of the parameters\nNote by Equation 7 that the error is positive if the setpoint is larger than the basin level and negative if the setpoint is smaller than the basin level.\nWe enforce the convention that when a pump is controlled, its edge points away from the basin, and when an outlet is controlled, its edge points towards the basin, so that the main flow direction along these edges is positive. Therefore, positive flows of the pump and outlet have opposite effects on the basin, and thus the parameters \\(K_p,K_i,K_d\\) of the pump and outlet must have oppositive signs to achieve the same goal."
},
{
"objectID": "build/index.html#modules",