diff --git a/.github/agents/critic.agent.md b/.github/agents/critic.agent.md new file mode 100644 index 00000000..fe6c6c5b --- /dev/null +++ b/.github/agents/critic.agent.md @@ -0,0 +1,21 @@ +--- +name: critic +description: Examines a Plugboard model and evaluates its suitability for a given task. +argument-hint: A description of the task the model is intended to solve, along with the model's design and implementation (e.g. a link to the code or documentation). +--- + +You are an independent critic responsible for evaluating the design of a Plugboard model. You will be given a description of the task the model is intended to solve, along with the model's design and implementation. Your job is to critically analyze the model and provide feedback on its strengths, weaknesses, and potential areas for improvement. + +## Your role +Examine the code you are provided. Think carefully and criticise the logic behind the model. Answer the following questions: +- Is it fit for purpose? +- What are its limitations? +- Does it use the plugboard framework appropriately? +- Is the code well-structured, readable and maintainable? Will it be clear and understandable to developers and data scientists who may need to work with it in the future? +- Has the model been implemented correctly? Review any assumptions and simplifications made in the model design, and evaluate whether they are reasonable given the task at hand. +- Examine any model output. Does it behave as expected? If not, can you identify why and suggest improvements to the model design or implementation that could help to address these issues? + +## Boundaries +- **Always** provide clear and constructive feedback that can help to improve the model. +- **Always** be specific in your critiques, citing particular aspects of the model's design or implementation. +- **If required** you can run the code and examine results to understand how it works. diff --git a/.github/agents/examples.agent.md b/.github/agents/examples.agent.md index 541334af..46c3f0e6 100644 --- a/.github/agents/examples.agent.md +++ b/.github/agents/examples.agent.md @@ -2,7 +2,7 @@ name: examples description: Develops example Plugboard models to demonstrate the capabilities of the framework argument-hint: A description of the example to generate, along with any specific requirements, ideas about structure, or constraints. -agents: ['researcher', 'docs', 'lint'] +agents: ['researcher', 'docs', 'lint', 'critic'] --- You are responsible for building high quality tutorials and demo examples for the Plugboard framework. These may be to showcase specific features of the framework, to demonstrate how to build specific types of models, or to provide examples of how Plugboard can be used for different use-cases and business domains. @@ -17,7 +17,6 @@ You are responsible for building high quality tutorials and demo examples for th - Prefer Jupyter notebooks for demo examples, as these allow for a mix of code, documentation and visualizations that can help to illustrate the concepts being demonstrated. - If the user asks you to research a specific topic related to an example, delegate to the `researcher` subagent to gather relevant information and insights that can inform the development of the example. - ## Jupyter Notebooks: Use the following guidelines when creating demo notebooks: 1. **Structure** @@ -34,10 +33,11 @@ Use the following guidelines when creating demo notebooks: - Include error handling 3. **Output** - Clear cell output before committing - - Generate plots where helpful + - Generate plots where helpful, prefer plotly for interactivity - Provide interpretation of results ## Boundaries: -- **Always** run the lint subagent on any code you write to ensure it adheres to the project's coding standards and is fully type-annotated. +- **Always** use the `critic` subagent to review your code and documentation, and make improvements based on its feedback to ensure that the examples are of high quality and behave as intended. +- **Always** run the `lint` subagent on any code you write to ensure it adheres to the project's coding standards and is fully type-annotated. - **Always** update the `mkdocs.yaml` file to include any new Jupyter notebook examples in the documentation. - **Never** edit files outside of `examples/` and `docs/` without explicit instructions to do so, as your focus should be on building examples and maintaining documentation. diff --git a/.github/agents/researcher.agent.md b/.github/agents/researcher.agent.md index 7ff3cb54..c28ee31f 100644 --- a/.github/agents/researcher.agent.md +++ b/.github/agents/researcher.agent.md @@ -2,7 +2,7 @@ name: researcher description: Researches specific topics on the internet and gathers relevant information. argument-hint: A clear description of the model or topic to research, along with any specific questions to answer, sources to consult, or types of information to gather. -tools: ['vscode', 'read', 'agent', 'search', 'web', 'todo'] +tools: ['vscode', 'read', 'agent', 'search', 'web', 'todo', 'notebooklm/*'] --- You are a subject-matter expert researcher responsible for gathering information on specific topics related to the Plugboard project. Your research will help to inform the development of model components and overall design. diff --git a/examples/demos/mining-minerals/001-crusher-circuit/.gitignore b/examples/demos/mining-minerals/001_crusher_circuit/.gitignore similarity index 100% rename from examples/demos/mining-minerals/001-crusher-circuit/.gitignore rename to examples/demos/mining-minerals/001_crusher_circuit/.gitignore diff --git a/examples/demos/mining-minerals/001-crusher-circuit/crusher_circuit.ipynb b/examples/demos/mining-minerals/001_crusher_circuit/crusher_circuit.ipynb similarity index 100% rename from examples/demos/mining-minerals/001-crusher-circuit/crusher_circuit.ipynb rename to examples/demos/mining-minerals/001_crusher_circuit/crusher_circuit.ipynb diff --git a/examples/demos/mining-minerals/001-crusher-circuit/crusher_circuit_components.py b/examples/demos/mining-minerals/001_crusher_circuit/crusher_circuit_components.py similarity index 100% rename from examples/demos/mining-minerals/001-crusher-circuit/crusher_circuit_components.py rename to examples/demos/mining-minerals/001_crusher_circuit/crusher_circuit_components.py diff --git a/examples/demos/physics-models/001-hot-water-tank/.gitignore b/examples/demos/physics-models/001_hot_water_tank/.gitignore similarity index 100% rename from examples/demos/physics-models/001-hot-water-tank/.gitignore rename to examples/demos/physics-models/001_hot_water_tank/.gitignore diff --git a/examples/demos/physics-models/001-hot-water-tank/hot-water-tank.ipynb b/examples/demos/physics-models/001_hot_water_tank/hot-water-tank.ipynb similarity index 100% rename from examples/demos/physics-models/001-hot-water-tank/hot-water-tank.ipynb rename to examples/demos/physics-models/001_hot_water_tank/hot-water-tank.ipynb diff --git a/examples/demos/physics-models/002_quadruple_tank/.gitignore b/examples/demos/physics-models/002_quadruple_tank/.gitignore new file mode 100644 index 00000000..d4621355 --- /dev/null +++ b/examples/demos/physics-models/002_quadruple_tank/.gitignore @@ -0,0 +1,2 @@ +*.csv +.ipynb_checkpoints/ diff --git a/examples/demos/physics-models/002_quadruple_tank/quadruple-tank.ipynb b/examples/demos/physics-models/002_quadruple_tank/quadruple-tank.ipynb new file mode 100644 index 00000000..5d59f774 --- /dev/null +++ b/examples/demos/physics-models/002_quadruple_tank/quadruple-tank.ipynb @@ -0,0 +1,1130 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Quadruple Tank Process\n", + "\n", + "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/plugboard-dev/plugboard/blob/main/examples/demos/physics-models/002-quadruple-tank/quadruple-tank.ipynb)\n", + "\n", + "[![Open in GitHub Codespaces](https://github.com/codespaces/badge.svg)](https://codespaces.new/plugboard-dev/plugboard)\n", + "\n", + "This notebook implements the **Quadruple Tank Process** — a classic MIMO (multi-input multi-output) benchmark problem from control theory, originally described by [Johansson (2000)](https://ieeexplore.ieee.org/document/845876).\n", + "\n", + "The system consists of four interconnected water tanks arranged in two columns:\n", + "- **Lower tanks** (Tank 1 & Tank 2): the controlled outputs\n", + "- **Upper tanks** (Tank 3 & Tank 4): fed by pumps, drain into the lower tanks in a cross-coupled fashion\n", + "- **Two pumps** with three-way valves split flow between upper and lower tanks\n", + "\n", + "By adjusting the valve split ratios ($\\gamma_1$, $\\gamma_2$), the system can exhibit either **minimum-phase** or **non-minimum-phase** behavior — a key property that fundamentally affects controllability.\n", + "\n", + "We'll build the nonlinear process model in Plugboard, add decentralized PI controllers, run closed-loop simulations, and visualize the results." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## Install and import dependencies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "# Install plugboard and dependencies for Google Colab\n", + "!pip install -q plugboard plotly" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "import typing as _t\n", + "\n", + "import pandas as pd\n", + "import plotly.graph_objects as go\n", + "from plotly.subplots import make_subplots\n", + "\n", + "from plugboard.component import Component\n", + "from plugboard.component import IOController as IO\n", + "from plugboard.connector import AsyncioConnector\n", + "from plugboard.library import FileWriter\n", + "from plugboard.process import LocalProcess\n", + "from plugboard.schemas import ComponentArgsDict, ConnectorSpec" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## Define quadruple tank process parameters\n", + "\n", + "We use the nominal parameter values from Johansson (2000). The key parameters are:\n", + "\n", + "| Parameter | Description | Value | Unit |\n", + "|-----------|-------------|-------|------|\n", + "| $A_1, A_3$ | Cross-section of Tanks 1, 3 | 28 | cm² |\n", + "| $A_2, A_4$ | Cross-section of Tanks 2, 4 | 32 | cm² |\n", + "| $a_1, a_3$ | Outlet hole area, Tanks 1, 3 | 0.071 | cm² |\n", + "| $a_2, a_4$ | Outlet hole area, Tanks 2, 4 | 0.057 | cm² |\n", + "| $k_1$ | Pump 1 flow constant | 3.33 | cm³/Vs |\n", + "| $k_2$ | Pump 2 flow constant | 3.35 | cm³/Vs |\n", + "| $g$ | Gravity | 981 | cm/s² |\n", + "| $\\gamma_1, \\gamma_2$ | Valve split ratios | varies | — |\n", + "\n", + "When $\\gamma_1 + \\gamma_2 > 1$, the system is **minimum phase** (easier to control).\n", + "When $\\gamma_1 + \\gamma_2 < 1$, the system is **non-minimum phase** (harder to control, exhibits inverse response)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "# Physical parameters (Johansson 2000)\n", + "A1 = 28.0 # cm² - cross-section of Tank 1\n", + "A2 = 32.0 # cm² - cross-section of Tank 2\n", + "A3 = 28.0 # cm² - cross-section of Tank 3\n", + "A4 = 32.0 # cm² - cross-section of Tank 4\n", + "\n", + "a1 = 0.071 # cm² - outlet hole area, Tank 1\n", + "a2 = 0.057 # cm² - outlet hole area, Tank 2\n", + "a3 = 0.071 # cm² - outlet hole area, Tank 3\n", + "a4 = 0.057 # cm² - outlet hole area, Tank 4\n", + "\n", + "g = 981.0 # cm/s² - gravitational acceleration\n", + "\n", + "k1 = 3.33 # cm³/Vs - pump 1 flow constant\n", + "k2 = 3.35 # cm³/Vs - pump 2 flow constant\n", + "\n", + "# Minimum-phase valve settings (gamma1 + gamma2 > 1)\n", + "GAMMA1_MP = 0.70\n", + "GAMMA2_MP = 0.60\n", + "\n", + "# Non-minimum-phase valve settings (gamma1 + gamma2 < 1)\n", + "GAMMA1_NMP = 0.43\n", + "GAMMA2_NMP = 0.34\n", + "\n", + "# Nominal operating point pump voltages\n", + "V1_0 = 3.0 # V\n", + "V2_0 = 3.0 # V\n", + "\n", + "# Steady-state levels for minimum-phase case (approximate, from Johansson 2000 Table I)\n", + "H1_0 = 12.4 # cm\n", + "H2_0 = 12.7 # cm\n", + "H3_0 = 1.8 # cm\n", + "H4_0 = 1.4 # cm\n", + "\n", + "# Steady-state levels for non-minimum-phase case.\n", + "# Upper tank steady state from Torricelli balance: h_i = ((1-γ)*k*V / (a*√(2g)))²\n", + "# Lower tank steady state includes drainage from the upper tank.\n", + "H3_NMP = ((1.0 - GAMMA2_NMP) * k2 * V2_0 / (a3 * math.sqrt(2.0 * g))) ** 2 # ≈ 4.45 cm\n", + "H4_NMP = ((1.0 - GAMMA1_NMP) * k1 * V1_0 / (a4 * math.sqrt(2.0 * g))) ** 2 # ≈ 5.09 cm\n", + "H1_NMP = (\n", + " (GAMMA1_NMP * k1 * V1_0 + a3 * math.sqrt(2.0 * g * H3_NMP)) / (a1 * math.sqrt(2.0 * g))\n", + ") ** 2\n", + "H2_NMP = (\n", + " (GAMMA2_NMP * k2 * V2_0 + a4 * math.sqrt(2.0 * g * H4_NMP)) / (a2 * math.sqrt(2.0 * g))\n", + ") ** 2\n", + "\n", + "# Simulation parameters\n", + "DT = 1.0 # s - time step\n", + "T_TOTAL = 600 # s - total simulation time\n", + "N_STEPS = int(T_TOTAL / DT)" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "## Build the tank components\n", + "\n", + "Each tank is modelled as an individual component, making the physical connections between tanks explicit in the process wiring. The nonlinear dynamics are governed by mass balance and Torricelli's law.\n", + "\n", + "**Lower tanks** (Tank 1, Tank 2) receive direct pump flow *and* drainage from the upper tanks:\n", + "\n", + "$$\\frac{dh_1}{dt} = -\\frac{a_1}{A_1}\\sqrt{2g\\,h_1} + \\frac{a_3}{A_1}\\sqrt{2g\\,h_3} + \\frac{\\gamma_1 k_1}{A_1}\\,v_1$$\n", + "\n", + "**Upper tanks** (Tank 3, Tank 4) receive only the cross-coupled pump flow:\n", + "\n", + "$$\\frac{dh_3}{dt} = -\\frac{a_3}{A_3}\\sqrt{2g\\,h_3} + \\frac{(1-\\gamma_2) k_2}{A_3}\\,v_2$$\n", + "\n", + "We use a simple Euler integration step with time step `dt`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "class LowerTank(Component):\n", + " \"\"\"Lower tank: receives direct pump flow and drainage from an upper tank.\"\"\"\n", + "\n", + " io = IO(inputs=[\"prev_h\", \"upper_h\", \"pump_voltage\"], outputs=[\"h\"])\n", + "\n", + " def __init__(\n", + " self,\n", + " tank_area: float = A1,\n", + " outlet_area: float = a1,\n", + " upper_outlet_area: float = a3,\n", + " pump_gain: float = k1,\n", + " gamma: float = GAMMA1_MP,\n", + " gravity: float = g,\n", + " dt: float = DT,\n", + " **kwargs: _t.Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._A = tank_area\n", + " self._a = outlet_area\n", + " self._a_upper = upper_outlet_area\n", + " self._k = pump_gain\n", + " self._gamma = gamma\n", + " self._dt = dt\n", + " self._sqrt_2g = math.sqrt(2.0 * gravity)\n", + "\n", + " async def step(self) -> None:\n", + " outflow = self._a * self._sqrt_2g * math.sqrt(max(self.prev_h, 0.0))\n", + " inflow_upper = self._a_upper * self._sqrt_2g * math.sqrt(max(self.upper_h, 0.0))\n", + " inflow_pump = self._gamma * self._k * self.pump_voltage\n", + "\n", + " dh = (-outflow + inflow_upper + inflow_pump) / self._A\n", + " self.h = max(self.prev_h + dh * self._dt, 0.0)\n", + "\n", + "\n", + "class UpperTank(Component):\n", + " \"\"\"Upper tank: receives cross-coupled pump flow only.\"\"\"\n", + "\n", + " io = IO(inputs=[\"prev_h\", \"pump_voltage\"], outputs=[\"h\"])\n", + "\n", + " def __init__(\n", + " self,\n", + " tank_area: float = A3,\n", + " outlet_area: float = a3,\n", + " pump_gain: float = k2,\n", + " gamma_complement: float = 1.0 - GAMMA2_MP,\n", + " gravity: float = g,\n", + " dt: float = DT,\n", + " **kwargs: _t.Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._A = tank_area\n", + " self._a = outlet_area\n", + " self._k = pump_gain\n", + " self._gamma_c = gamma_complement\n", + " self._dt = dt\n", + " self._sqrt_2g = math.sqrt(2.0 * gravity)\n", + "\n", + " async def step(self) -> None:\n", + " outflow = self._a * self._sqrt_2g * math.sqrt(max(self.prev_h, 0.0))\n", + " inflow_pump = self._gamma_c * self._k * self.pump_voltage\n", + "\n", + " dh = (-outflow + inflow_pump) / self._A\n", + " self.h = max(self.prev_h + dh * self._dt, 0.0)" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "## Build the PI Controller component\n", + "\n", + "We use **decentralized PI control** — two independent loops:\n", + "- Loop 1: $v_1$ controls $h_1$\n", + "- Loop 2: $v_2$ controls $h_2$\n", + "\n", + "Each controller computes:\n", + "$$u(t) = v_0 + K_p \\cdot e(t) + K_i \\int_0^t e(\\tau)\\,d\\tau$$\n", + "\n", + "where $v_0$ is the steady-state pump voltage (feedforward), and the PI terms provide feedback correction. Anti-windup clamping limits control output between 0V and 10V." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "class PIController(Component):\n", + " \"\"\"Decentralized PI controller with anti-windup for two tank levels.\"\"\"\n", + "\n", + " io = IO(\n", + " inputs=[\"h1\", \"h2\", \"sp_h1\", \"sp_h2\"],\n", + " outputs=[\"v1\", \"v2\"],\n", + " )\n", + "\n", + " def __init__(\n", + " self,\n", + " kp: tuple[float, float] = (3.0, 2.7),\n", + " ki: tuple[float, float] = (0.1, 0.068),\n", + " v_min: float = 0.0,\n", + " v_max: float = 10.0,\n", + " v0: tuple[float, float] = (V1_0, V2_0),\n", + " dt: float = DT,\n", + " **kwargs: _t.Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._kp = kp\n", + " self._ki = ki\n", + " self._v_min = v_min\n", + " self._v_max = v_max\n", + " self._v0 = v0\n", + " self._dt = dt\n", + " self._integral = [0.0, 0.0]\n", + "\n", + " async def step(self) -> None:\n", + " errors = [self.sp_h1 - self.h1, self.sp_h2 - self.h2]\n", + "\n", + " outputs = []\n", + " for i, (kp, ki, v0, e) in enumerate(zip(self._kp, self._ki, self._v0, errors)):\n", + " u_raw = v0 + kp * e + ki * self._integral[i]\n", + " u_clamped = max(self._v_min, min(self._v_max, u_raw))\n", + " # Anti-windup: only integrate when the output is not saturated\n", + " if u_raw == u_clamped:\n", + " self._integral[i] += e * self._dt\n", + " outputs.append(u_clamped)\n", + "\n", + " # PI output = feedforward (steady-state voltage) + feedback\n", + " self.v1, self.v2 = outputs[0], outputs[1]" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "## Build the Setpoint Generator component\n", + "\n", + "The setpoint generator produces step changes:\n", + "- $t=0$: both setpoints at steady-state levels ($h_1^0 = 12.4$ cm, $h_2^0 = 12.7$ cm)\n", + "- $t=100$ s: step increase in $h_1$ setpoint by +2 cm\n", + "- $t=300$ s: step increase in $h_2$ setpoint by +2 cm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "class SetpointGenerator(Component):\n", + " \"\"\"Generates step changes in setpoints for h1 and h2.\"\"\"\n", + "\n", + " io = IO(outputs=[\"sp_h1\", \"sp_h2\"])\n", + "\n", + " def __init__(\n", + " self,\n", + " h1_base: float = H1_0,\n", + " h2_base: float = H2_0,\n", + " h1_step: float = 2.0,\n", + " h2_step: float = 2.0,\n", + " h1_step_time: int = 100,\n", + " h2_step_time: int = 300,\n", + " n_steps: int = N_STEPS,\n", + " **kwargs: _t.Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._h1_base = h1_base\n", + " self._h2_base = h2_base\n", + " self._h1_step = h1_step\n", + " self._h2_step = h2_step\n", + " self._h1_step_time = h1_step_time\n", + " self._h2_step_time = h2_step_time\n", + " self._n_steps = n_steps\n", + " self._iteration = 0\n", + "\n", + " async def step(self) -> None:\n", + " t = self._iteration\n", + "\n", + " self.sp_h1 = self._h1_base + (self._h1_step if t >= self._h1_step_time else 0.0)\n", + " self.sp_h2 = self._h2_base + (self._h2_step if t >= self._h2_step_time else 0.0)\n", + "\n", + " self._iteration += 1\n", + " if self._iteration >= self._n_steps:\n", + " await self.io.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "def build_process(\n", + " gamma1: float,\n", + " gamma2: float,\n", + " kp: tuple[float, float],\n", + " ki: tuple[float, float],\n", + " csv_path: str,\n", + " h1_init: float,\n", + " h2_init: float,\n", + " h3_init: float,\n", + " h4_init: float,\n", + ") -> LocalProcess:\n", + " \"\"\"Factory for a closed-loop quadruple-tank process.\n", + "\n", + " All connector wiring is identical across scenarios; only the valve ratios,\n", + " controller gains, and initial conditions differ.\n", + " \"\"\"\n", + "\n", + " def _connect(src: str, tgt: str) -> AsyncioConnector:\n", + " return AsyncioConnector(spec=ConnectorSpec(source=src, target=tgt))\n", + "\n", + " comps = [\n", + " SetpointGenerator(name=\"setpoints\"),\n", + " PIController(\n", + " name=\"controller\", kp=kp, ki=ki, initial_values={\"h1\": [h1_init], \"h2\": [h2_init]}\n", + " ),\n", + " LowerTank(\n", + " name=\"tank1\",\n", + " tank_area=A1,\n", + " outlet_area=a1,\n", + " upper_outlet_area=a3,\n", + " pump_gain=k1,\n", + " gamma=gamma1,\n", + " initial_values={\"prev_h\": [h1_init]},\n", + " ),\n", + " LowerTank(\n", + " name=\"tank2\",\n", + " tank_area=A2,\n", + " outlet_area=a2,\n", + " upper_outlet_area=a4,\n", + " pump_gain=k2,\n", + " gamma=gamma2,\n", + " initial_values={\"prev_h\": [h2_init]},\n", + " ),\n", + " UpperTank(\n", + " name=\"tank3\",\n", + " tank_area=A3,\n", + " outlet_area=a3,\n", + " pump_gain=k2,\n", + " gamma_complement=1.0 - gamma2,\n", + " initial_values={\"prev_h\": [h3_init]},\n", + " ),\n", + " UpperTank(\n", + " name=\"tank4\",\n", + " tank_area=A4,\n", + " outlet_area=a4,\n", + " pump_gain=k1,\n", + " gamma_complement=1.0 - gamma1,\n", + " initial_values={\"prev_h\": [h4_init]},\n", + " ),\n", + " FileWriter(\n", + " name=\"save\",\n", + " path=csv_path,\n", + " field_names=[\"h1\", \"h2\", \"h3\", \"h4\", \"v1\", \"v2\", \"sp_h1\", \"sp_h2\"],\n", + " ),\n", + " ]\n", + "\n", + " conns = [\n", + " # Setpoints → Controller\n", + " _connect(\"setpoints.sp_h1\", \"controller.sp_h1\"),\n", + " _connect(\"setpoints.sp_h2\", \"controller.sp_h2\"),\n", + " # Controller → Pumps → Tanks (direct and cross-coupled flow)\n", + " _connect(\"controller.v1\", \"tank1.pump_voltage\"), # Pump 1 → Tank 1 (lower left)\n", + " _connect(\"controller.v2\", \"tank2.pump_voltage\"), # Pump 2 → Tank 2 (lower right)\n", + " _connect(\n", + " \"controller.v2\", \"tank3.pump_voltage\"\n", + " ), # Pump 2 → Tank 3 (upper left, cross-coupled)\n", + " _connect(\n", + " \"controller.v1\", \"tank4.pump_voltage\"\n", + " ), # Pump 1 → Tank 4 (upper right, cross-coupled)\n", + " # Upper tanks drain into lower tanks\n", + " _connect(\"tank3.h\", \"tank1.upper_h\"),\n", + " _connect(\"tank4.h\", \"tank2.upper_h\"),\n", + " # Tank feedback → Controller\n", + " _connect(\"tank1.h\", \"controller.h1\"),\n", + " _connect(\"tank2.h\", \"controller.h2\"),\n", + " # Self-feedback for Euler state (h at step t becomes prev_h at step t+1)\n", + " _connect(\"tank1.h\", \"tank1.prev_h\"),\n", + " _connect(\"tank2.h\", \"tank2.prev_h\"),\n", + " _connect(\"tank3.h\", \"tank3.prev_h\"),\n", + " _connect(\"tank4.h\", \"tank4.prev_h\"),\n", + " # Save all signals\n", + " _connect(\"tank1.h\", \"save.h1\"),\n", + " _connect(\"tank2.h\", \"save.h2\"),\n", + " _connect(\"tank3.h\", \"save.h3\"),\n", + " _connect(\"tank4.h\", \"save.h4\"),\n", + " _connect(\"controller.v1\", \"save.v1\"),\n", + " _connect(\"controller.v2\", \"save.v2\"),\n", + " _connect(\"setpoints.sp_h1\", \"save.sp_h1\"),\n", + " _connect(\"setpoints.sp_h2\", \"save.sp_h2\"),\n", + " ]\n", + "\n", + " return LocalProcess(components=comps, connectors=conns)" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "## Assemble the process\n", + "\n", + "Each tank is a separate component, so the physical coupling between them is explicit in the connector wiring.\n", + "\n", + "Cross-coupling: Tank 3 (upper left) drains into Tank 1 (lower left), Tank 4 (upper right) drains into Tank 2 (lower right). Notice how Pump 1 feeds Tank 1 directly *and* Tank 4 via the cross-coupled valve.\n", + "\n", + "**State propagation via self-feedback**: Plugboard connectors deliver values with a one-step lag. We exploit this for Euler integration: `tank1.h → tank1.prev_h` routes each step's computed height back as the next step's starting value. The `initial_values` parameter seeds `prev_h` for the very first iteration." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "# Minimum-phase simulation using nominal Johansson (2000) parameter values\n", + "process = build_process(\n", + " gamma1=GAMMA1_MP,\n", + " gamma2=GAMMA2_MP,\n", + " kp=(3.0, 2.7),\n", + " ki=(0.1, 0.068),\n", + " csv_path=\"quadruple_tank_mp.csv\",\n", + " h1_init=H1_0,\n", + " h2_init=H2_0,\n", + " h3_init=H3_0,\n", + " h4_init=H4_0,\n", + ")\n", + "print(f\"{len(process.components)} components, {len(process.connectors)} connectors\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "from IPython.display import Image\n", + "from plugboard.diagram import MermaidDiagram\n", + "\n", + "Image(url=MermaidDiagram.from_process(process).url)" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "## Run the simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "async with process:\n", + " await process.run()\n", + "\n", + "df_mp = pd.read_csv(\"quadruple_tank_mp.csv\")\n", + "df_mp[\"time\"] = df_mp.index * DT\n", + "df_mp.head()" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "## Plot tank levels over time\n", + "\n", + "The top subplot shows the controlled outputs ($h_1$, $h_2$) and their setpoints.\n", + "The bottom subplot shows the upper tank levels ($h_3$, $h_4$)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "fig = make_subplots(\n", + " rows=2,\n", + " cols=1,\n", + " shared_xaxes=True,\n", + " subplot_titles=(\"Lower tank levels (controlled)\", \"Upper tank levels\"),\n", + " vertical_spacing=0.1,\n", + ")\n", + "\n", + "# Lower tank levels and setpoints\n", + "fig.add_trace(\n", + " go.Scatter(x=df_mp[\"time\"], y=df_mp[\"h1\"], name=\"h₁\", line=dict(color=\"blue\")), row=1, col=1\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(x=df_mp[\"time\"], y=df_mp[\"h2\"], name=\"h₂\", line=dict(color=\"red\")), row=1, col=1\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_mp[\"time\"], y=df_mp[\"sp_h1\"], name=\"SP h₁\", line=dict(color=\"blue\", dash=\"dash\")\n", + " ),\n", + " row=1,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_mp[\"time\"], y=df_mp[\"sp_h2\"], name=\"SP h₂\", line=dict(color=\"red\", dash=\"dash\")\n", + " ),\n", + " row=1,\n", + " col=1,\n", + ")\n", + "\n", + "# Upper tank levels\n", + "fig.add_trace(\n", + " go.Scatter(x=df_mp[\"time\"], y=df_mp[\"h3\"], name=\"h₃\", line=dict(color=\"green\")), row=2, col=1\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(x=df_mp[\"time\"], y=df_mp[\"h4\"], name=\"h₄\", line=dict(color=\"orange\")), row=2, col=1\n", + ")\n", + "\n", + "fig.update_layout(height=600, title_text=\"Quadruple Tank Process — Minimum Phase (γ₁=0.7, γ₂=0.6)\")\n", + "fig.update_xaxes(title_text=\"Time (s)\", row=2, col=1)\n", + "fig.update_yaxes(title_text=\"Level (cm)\", row=1, col=1)\n", + "fig.update_yaxes(title_text=\"Level (cm)\", row=2, col=1)\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "## Plot control inputs\n", + "\n", + "This shows how the pump voltages $v_1$ and $v_2$ respond to setpoint changes. The dashed lines indicate the actuator saturation limits (0V and 10V)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "fig = go.Figure()\n", + "fig.add_trace(go.Scatter(x=df_mp[\"time\"], y=df_mp[\"v1\"], name=\"v₁\", line=dict(color=\"blue\")))\n", + "fig.add_trace(go.Scatter(x=df_mp[\"time\"], y=df_mp[\"v2\"], name=\"v₂\", line=dict(color=\"red\")))\n", + "fig.add_hline(y=0.0, line_dash=\"dash\", line_color=\"gray\", annotation_text=\"Min (0V)\")\n", + "fig.add_hline(y=10.0, line_dash=\"dash\", line_color=\"gray\", annotation_text=\"Max (10V)\")\n", + "fig.update_layout(\n", + " title=\"Pump Voltages — Minimum Phase\",\n", + " xaxis_title=\"Time (s)\",\n", + " yaxis_title=\"Voltage (V)\",\n", + " height=400,\n", + ")\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "## Setpoint tracking performance\n", + "\n", + "We compute the tracking error for each loop and report key performance metrics." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23", + "metadata": {}, + "outputs": [], + "source": [ + "df_mp[\"e1\"] = df_mp[\"sp_h1\"] - df_mp[\"h1\"]\n", + "df_mp[\"e2\"] = df_mp[\"sp_h2\"] - df_mp[\"h2\"]\n", + "\n", + "fig = go.Figure()\n", + "fig.add_trace(go.Scatter(x=df_mp[\"time\"], y=df_mp[\"e1\"], name=\"Error h₁\", line=dict(color=\"blue\")))\n", + "fig.add_trace(go.Scatter(x=df_mp[\"time\"], y=df_mp[\"e2\"], name=\"Error h₂\", line=dict(color=\"red\")))\n", + "fig.add_hline(y=0.0, line_dash=\"dash\", line_color=\"gray\")\n", + "fig.update_layout(\n", + " title=\"Tracking Error — Minimum Phase\",\n", + " xaxis_title=\"Time (s)\",\n", + " yaxis_title=\"Error (cm)\",\n", + " height=400,\n", + ")\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24", + "metadata": {}, + "outputs": [], + "source": [ + "def compute_metrics(\n", + " df: pd.DataFrame, sp_col: str, val_col: str, step_time: float\n", + ") -> dict[str, float]:\n", + " \"\"\"Compute step response metrics for a single control loop.\n", + "\n", + " Parameters\n", + " ----------\n", + " step_time:\n", + " Time (in seconds, matching ``df[\"time\"]``) at which the step change occurs.\n", + " \"\"\"\n", + " mask = df[\"time\"] >= step_time\n", + " subset = df.loc[mask].copy()\n", + "\n", + " if subset.empty:\n", + " return {\n", + " \"rise_time\": float(\"nan\"),\n", + " \"settling_time_from_step\": float(\"nan\"),\n", + " \"overshoot_pct\": float(\"nan\"),\n", + " \"ss_error\": float(\"nan\"),\n", + " }\n", + "\n", + " try:\n", + " sp_final = float(subset[sp_col].iloc[-1])\n", + " except (KeyError, IndexError, ValueError):\n", + " return {\n", + " \"rise_time\": float(\"nan\"),\n", + " \"settling_time_from_step\": float(\"nan\"),\n", + " \"overshoot_pct\": float(\"nan\"),\n", + " \"ss_error\": float(\"nan\"),\n", + " }\n", + "\n", + " try:\n", + " sp_initial = float(df.loc[df[\"time\"] < step_time, sp_col].iloc[-1])\n", + " except (IndexError, KeyError):\n", + " sp_initial = float(subset[sp_col].iloc[0])\n", + "\n", + " step_size = sp_final - sp_initial\n", + "\n", + " if abs(step_size) < 1e-6:\n", + " return {\n", + " \"rise_time\": float(\"nan\"),\n", + " \"settling_time_from_step\": float(\"nan\"),\n", + " \"overshoot_pct\": float(\"nan\"),\n", + " \"ss_error\": float(\"nan\"),\n", + " }\n", + "\n", + " direction = 1 if step_size > 0 else -1\n", + "\n", + " val_10 = sp_initial + 0.1 * step_size\n", + " val_90 = sp_initial + 0.9 * step_size\n", + "\n", + " if direction > 0:\n", + " cond_10 = subset[val_col] >= val_10\n", + " cond_90 = subset[val_col] >= val_90\n", + " else:\n", + " cond_10 = subset[val_col] <= val_10\n", + " cond_90 = subset[val_col] <= val_90\n", + "\n", + " t_10 = subset.loc[cond_10, \"time\"].iloc[0] if cond_10.any() else float(\"nan\")\n", + " t_90 = subset.loc[cond_90, \"time\"].iloc[0] if cond_90.any() else float(\"nan\")\n", + "\n", + " rise_time = t_90 - t_10 if not (math.isnan(t_10) or math.isnan(t_90)) else float(\"nan\")\n", + "\n", + " # Settling time: last time error exceeds 2% of step size, measured from the step change\n", + " error = abs(subset[val_col] - sp_final)\n", + " threshold = 0.02 * abs(step_size)\n", + " settled = error <= threshold\n", + " if settled.all():\n", + " settling_time = 0.0\n", + " elif settled.any():\n", + " # last time signal is outside tolerance\n", + " last_out_of_bound = subset.loc[~settled, \"time\"].iloc[-1]\n", + " settling_time = last_out_of_bound - step_time\n", + " else:\n", + " settling_time = float(\"nan\")\n", + "\n", + " # Overshoot: consider direction of step\n", + " peak = subset[val_col].max() if direction > 0 else subset[val_col].min()\n", + " if direction > 0:\n", + " overshoot_pct = max(0.0, (peak - sp_final) / abs(step_size) * 100.0)\n", + " else:\n", + " overshoot_pct = max(0.0, (sp_final - peak) / abs(step_size) * 100.0)\n", + "\n", + " # Steady-state error (average over last 10% of data)\n", + " tail = subset.tail(max(1, len(subset) // 10))\n", + " ss_error = abs(tail[sp_col].mean() - tail[val_col].mean())\n", + "\n", + " return {\n", + " \"rise_time\": round(rise_time, 1) if not math.isnan(rise_time) else float(\"nan\"),\n", + " \"settling_time_from_step\": round(settling_time, 1)\n", + " if not math.isnan(settling_time)\n", + " else float(\"nan\"),\n", + " \"overshoot_pct\": round(overshoot_pct, 1),\n", + " \"ss_error\": round(ss_error, 3),\n", + " }\n", + "\n", + "\n", + "metrics_h1 = compute_metrics(df_mp, \"sp_h1\", \"h1\", 100.0)\n", + "metrics_h2 = compute_metrics(df_mp, \"sp_h2\", \"h2\", 300.0)\n", + "\n", + "pd.DataFrame({\"h₁ loop\": metrics_h1, \"h₂ loop\": metrics_h2}).rename(\n", + " index={\n", + " \"rise_time\": \"Rise time (s)\",\n", + " \"settling_time_from_step\": \"Settling time from step (s)\",\n", + " \"overshoot_pct\": \"Overshoot (%)\",\n", + " \"ss_error\": \"Steady-state error (cm)\",\n", + " },\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "25", + "metadata": {}, + "source": [ + "## Minimum-phase vs non-minimum-phase\n", + "\n", + "The key insight of the quadruple tank process is that valve parameters $\\gamma_1$ and $\\gamma_2$ determine the system's phase characteristics:\n", + "\n", + "- **Minimum phase** ($\\gamma_1 + \\gamma_2 > 1$): Most of the pump flow goes directly to the lower tanks. The system responds promptly to control inputs.\n", + "- **Non-minimum phase** ($\\gamma_1 + \\gamma_2 < 1$): Most of the flow goes to the upper tanks first, causing **inverse response** — the level initially moves in the wrong direction before correcting. This fundamentally limits achievable control performance.\n", + "\n", + "Let's re-run the simulation with non-minimum-phase parameters and compare." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26", + "metadata": {}, + "outputs": [], + "source": [ + "# Non-minimum-phase simulation — uses correct NMP steady-state initial conditions\n", + "# and more conservative PI gains (NMP systems require reduced bandwidth)\n", + "process_nmp = build_process(\n", + " gamma1=GAMMA1_NMP,\n", + " gamma2=GAMMA2_NMP,\n", + " kp=(0.5, 0.5),\n", + " ki=(0.005, 0.004),\n", + " csv_path=\"quadruple_tank_nmp.csv\",\n", + " h1_init=H1_NMP,\n", + " h2_init=H2_NMP,\n", + " h3_init=H3_NMP,\n", + " h4_init=H4_NMP,\n", + ")\n", + "\n", + "async with process_nmp:\n", + " await process_nmp.run()\n", + "\n", + "df_nmp = pd.read_csv(\"quadruple_tank_nmp.csv\")\n", + "df_nmp[\"time\"] = df_nmp.index * DT" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27", + "metadata": {}, + "outputs": [], + "source": [ + "fig = make_subplots(\n", + " rows=2,\n", + " cols=2,\n", + " shared_xaxes=True,\n", + " subplot_titles=(\n", + " \"Minimum phase (γ₁=0.7, γ₂=0.6)\",\n", + " \"Non-minimum phase (γ₁=0.43, γ₂=0.34)\",\n", + " \"Upper tanks (MP)\",\n", + " \"Upper tanks (NMP)\",\n", + " ),\n", + " vertical_spacing=0.12,\n", + " horizontal_spacing=0.08,\n", + ")\n", + "\n", + "# Minimum phase — lower tanks\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_mp[\"time\"], y=df_mp[\"h1\"], name=\"h₁ (MP)\", line=dict(color=\"blue\"), legendgroup=\"mp\"\n", + " ),\n", + " row=1,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_mp[\"time\"], y=df_mp[\"h2\"], name=\"h₂ (MP)\", line=dict(color=\"red\"), legendgroup=\"mp\"\n", + " ),\n", + " row=1,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_mp[\"time\"],\n", + " y=df_mp[\"sp_h1\"],\n", + " name=\"SP h₁\",\n", + " line=dict(color=\"blue\", dash=\"dash\"),\n", + " legendgroup=\"mp\",\n", + " showlegend=False,\n", + " ),\n", + " row=1,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_mp[\"time\"],\n", + " y=df_mp[\"sp_h2\"],\n", + " name=\"SP h₂\",\n", + " line=dict(color=\"red\", dash=\"dash\"),\n", + " legendgroup=\"mp\",\n", + " showlegend=False,\n", + " ),\n", + " row=1,\n", + " col=1,\n", + ")\n", + "\n", + "# Non-minimum phase — lower tanks\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_nmp[\"time\"],\n", + " y=df_nmp[\"h1\"],\n", + " name=\"h₁ (NMP)\",\n", + " line=dict(color=\"blue\", dash=\"dot\"),\n", + " legendgroup=\"nmp\",\n", + " ),\n", + " row=1,\n", + " col=2,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_nmp[\"time\"],\n", + " y=df_nmp[\"h2\"],\n", + " name=\"h₂ (NMP)\",\n", + " line=dict(color=\"red\", dash=\"dot\"),\n", + " legendgroup=\"nmp\",\n", + " ),\n", + " row=1,\n", + " col=2,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_nmp[\"time\"],\n", + " y=df_nmp[\"sp_h1\"],\n", + " name=\"SP h₁\",\n", + " line=dict(color=\"blue\", dash=\"dash\"),\n", + " legendgroup=\"nmp\",\n", + " showlegend=False,\n", + " ),\n", + " row=1,\n", + " col=2,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_nmp[\"time\"],\n", + " y=df_nmp[\"sp_h2\"],\n", + " name=\"SP h₂\",\n", + " line=dict(color=\"red\", dash=\"dash\"),\n", + " legendgroup=\"nmp\",\n", + " showlegend=False,\n", + " ),\n", + " row=1,\n", + " col=2,\n", + ")\n", + "\n", + "# Minimum phase — upper tanks\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_mp[\"time\"], y=df_mp[\"h3\"], name=\"h₃ (MP)\", line=dict(color=\"green\"), legendgroup=\"mp\"\n", + " ),\n", + " row=2,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_mp[\"time\"], y=df_mp[\"h4\"], name=\"h₄ (MP)\", line=dict(color=\"orange\"), legendgroup=\"mp\"\n", + " ),\n", + " row=2,\n", + " col=1,\n", + ")\n", + "\n", + "# Non-minimum phase — upper tanks\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_nmp[\"time\"],\n", + " y=df_nmp[\"h3\"],\n", + " name=\"h₃ (NMP)\",\n", + " line=dict(color=\"green\", dash=\"dot\"),\n", + " legendgroup=\"nmp\",\n", + " ),\n", + " row=2,\n", + " col=2,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_nmp[\"time\"],\n", + " y=df_nmp[\"h4\"],\n", + " name=\"h₄ (NMP)\",\n", + " line=dict(color=\"orange\", dash=\"dot\"),\n", + " legendgroup=\"nmp\",\n", + " ),\n", + " row=2,\n", + " col=2,\n", + ")\n", + "\n", + "fig.update_layout(height=700, title_text=\"Minimum Phase vs Non-Minimum Phase Comparison\")\n", + "fig.update_xaxes(title_text=\"Time (s)\", row=2, col=1)\n", + "fig.update_xaxes(title_text=\"Time (s)\", row=2, col=2)\n", + "for row in [1, 2]:\n", + " for col in [1, 2]:\n", + " fig.update_yaxes(title_text=\"Level (cm)\", row=row, col=col)\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28", + "metadata": {}, + "outputs": [], + "source": [ + "# Zoom in on the NMP h1 response to reveal the inverse‑response (wrong‑direction) transient.\n", + "# For a non-minimum-phase system the step in sp_h2 at t=300 s causes h1 to dip *down*\n", + "# before recovering — a hallmark characteristic of cross-coupled NMP dynamics.\n", + "ZOOM_START, ZOOM_END = 295, 360\n", + "\n", + "mask = (df_nmp[\"time\"] >= ZOOM_START) & (df_nmp[\"time\"] <= ZOOM_END)\n", + "zoom = df_nmp.loc[mask]\n", + "\n", + "fig_zoom = go.Figure()\n", + "fig_zoom.add_trace(\n", + " go.Scatter(\n", + " x=zoom[\"time\"],\n", + " y=zoom[\"h1\"],\n", + " name=\"h₁ (NMP)\",\n", + " line=dict(color=\"blue\"),\n", + " )\n", + ")\n", + "fig_zoom.add_trace(\n", + " go.Scatter(\n", + " x=zoom[\"time\"],\n", + " y=zoom[\"sp_h1\"],\n", + " name=\"SP h₁\",\n", + " line=dict(color=\"blue\", dash=\"dash\"),\n", + " )\n", + ")\n", + "# Mark the step change\n", + "fig_zoom.add_vline(\n", + " x=300,\n", + " line_dash=\"dot\",\n", + " line_color=\"gray\",\n", + " annotation_text=\"sp_h2 step\",\n", + " annotation_position=\"top right\",\n", + ")\n", + "# Annotate the inverse-response dip\n", + "inv_idx = zoom[\"h1\"].idxmin()\n", + "fig_zoom.add_annotation(\n", + " x=zoom.loc[inv_idx, \"time\"],\n", + " y=zoom.loc[inv_idx, \"h1\"],\n", + " text=\"Inverse response
(initial decrease)\",\n", + " showarrow=True,\n", + " arrowhead=2,\n", + " ax=40,\n", + " ay=-40,\n", + " font=dict(size=12),\n", + ")\n", + "fig_zoom.update_layout(\n", + " title=\"NMP inverse response — zoom around h₂ step change (t = 300 s)\",\n", + " xaxis_title=\"Time (s)\",\n", + " yaxis_title=\"Tank level (cm)\",\n", + " legend=dict(orientation=\"h\", y=-0.2),\n", + " template=\"plotly_white\",\n", + ")\n", + "fig_zoom.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29", + "metadata": {}, + "outputs": [], + "source": [ + "# Compare performance metrics\n", + "metrics_nmp_h1 = compute_metrics(df_nmp, \"sp_h1\", \"h1\", 100.0)\n", + "metrics_nmp_h2 = compute_metrics(df_nmp, \"sp_h2\", \"h2\", 300.0)\n", + "\n", + "pd.DataFrame(\n", + " {\n", + " \"h₁ (MP)\": metrics_h1,\n", + " \"h₁ (NMP)\": metrics_nmp_h1,\n", + " \"h₂ (MP)\": metrics_h2,\n", + " \"h₂ (NMP)\": metrics_nmp_h2,\n", + " }\n", + ").rename(\n", + " index={\n", + " \"rise_time\": \"Rise time (s)\",\n", + " \"settling_time_from_step\": \"Settling time from step (s)\",\n", + " \"overshoot_pct\": \"Overshoot (%)\",\n", + " \"ss_error\": \"Steady-state error (cm)\",\n", + " }\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "30", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "This notebook demonstrated how to model and simulate the **quadruple tank process** using Plugboard:\n", + "\n", + "1. **Individual tank components**: Each tank is a separate `LowerTank` or `UpperTank` component, making the physical connections between tanks explicit in the process wiring — you can directly see which pump feeds which tank, and how upper tanks drain into lower tanks.\n", + "2. **Controller**: A decentralized PI controller with anti-windup clamping regulates the lower tank levels.\n", + "3. **Setpoint generator**: Step changes in both setpoints showcase the MIMO coupling between loops.\n", + "4. **Phase behavior**: By changing the valve split ratios ($\\gamma_1$, $\\gamma_2$), we demonstrated the fundamental difference between minimum-phase and non-minimum-phase configurations — the latter exhibits slower response due to inverse response dynamics.\n", + "\n", + "### Key takeaways\n", + "- Modelling each tank as a separate component makes the cross-coupling visible: Pump 1 feeds Tank 1 directly *and* Tank 4 via the valve, while Tank 3 drains into Tank 1.\n", + "- The minimum-phase configuration responds faster because pump flow goes directly to the controlled (lower) tanks.\n", + "- The non-minimum-phase configuration requires more conservative tuning and responds much more slowly, because the flow must pass through the upper tanks first.\n", + "- This benchmark illustrates why MIMO systems with right-half-plane zeros pose fundamental performance limitations for feedback control.\n", + "\n", + "### References\n", + "- Johansson, K.H. (2000). \"The Quadruple-Tank Process: A Multivariable Laboratory Process with an Adjustable Zero.\" *IEEE Transactions on Control Systems Technology*, 8(3), pp. 456–465." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "plugboard (3.12.8)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/demos/process-control/001_bsm1_wastewater/.gitignore b/examples/demos/process-control/001_bsm1_wastewater/.gitignore new file mode 100644 index 00000000..a4b0eb75 --- /dev/null +++ b/examples/demos/process-control/001_bsm1_wastewater/.gitignore @@ -0,0 +1 @@ +bsm1_results.csv diff --git a/examples/demos/process-control/001_bsm1_wastewater/bsm1-wastewater.ipynb b/examples/demos/process-control/001_bsm1_wastewater/bsm1-wastewater.ipynb new file mode 100644 index 00000000..ef9355dc --- /dev/null +++ b/examples/demos/process-control/001_bsm1_wastewater/bsm1-wastewater.ipynb @@ -0,0 +1,1203 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Reduced BSM1-Inspired Wastewater Treatment Model\n", + "\n", + "\n", + "\n", + "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/plugboard-dev/plugboard/blob/main/examples/demos/process-control/001_bsm1_wastewater/bsm1-wastewater.ipynb)\n", + "\n", + "\n", + "\n", + "[![Open in GitHub Codespaces](https://github.com/codespaces/badge.svg)](https://codespaces.new/plugboard-dev/plugboard)\n", + "\n", + "\n", + "\n", + "This notebook implements a reduced BSM1-inspired activated sludge plant in Plugboard. BSM1, the Benchmark Simulation Model No. 1, is a standard wastewater-process benchmark built around ASM1 biological kinetics, pre-denitrification, and closed-loop control.\n", + "\n", + "\n", + "\n", + "The full benchmark is too large and stiff for a clear teaching notebook, so this example focuses on the parts that best demonstrate Plugboard:\n", + "\n", + "- one anoxic reactor and one aerobic reactor\n", + "\n", + "- an internal recycle loop for nitrate return\n", + "\n", + "- a fixed return sludge stream from an ideal settler\n", + "\n", + "- PI control of nitrate and dissolved oxygen\n", + "\n", + "- dynamic influent disturbances, file-backed simulation output, and visualization\n", + "\n", + "\n", + "\n", + "This is a control-topology tutorial, not a benchmark-faithful BSM1 reproduction. The notebook keeps the broad pre-denitrification story and the controller roles, but it also uses a reduced ASM1-inspired state vector and surrogate kinetics inside the reactors. The reported values are therefore illustrative tutorial outputs rather than benchmark-comparable KPIs, controller-tuning guidance, or publication-grade plant predictions.\n", + "\n", + "\n", + "\n", + "The component classes and process assembly are shown directly in the notebook so readers can inspect the Plugboard model structure cell by cell." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## Install and import dependencies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "# Install plugboard and plotting dependencies for Google Colab\n", + "%pip install -q plugboard plotly pandas" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "import math\n", + "import typing as _t\n", + "\n", + "import pandas as pd\n", + "import plotly.graph_objects as go\n", + "from plotly.subplots import make_subplots\n", + "from pydantic import BaseModel, ConfigDict\n", + "\n", + "from plugboard.component import Component\n", + "from plugboard.component import IOController as IO\n", + "from plugboard.connector import AsyncioConnector\n", + "from plugboard.diagram import MermaidDiagram, markdown_diagram\n", + "from plugboard.library import FileWriter\n", + "from plugboard.process import LocalProcess\n", + "from plugboard.schemas import ComponentArgsDict, ConnectorSpec" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## Which BSM1 ideas are preserved, reduced, or not included?\n", + "\n", + "The table below makes the scope of the model explicit so the reader does not confuse this demo with a benchmark-grade BSM1 implementation. It preserves the plant architecture and control roles at a high level, but not the full benchmark dynamics or KPI machinery.\n", + "\n", + "| Status | Aspect | Explanation |\n", + "| --- | --- | --- |\n", + "| Preserved | Two PI loops | DO is manipulated through KLa and nitrate is manipulated through internal recycle flow, so the controller roles still match the standard BSM1 control story. |\n", + "| Preserved | Fixed return sludge flow | The return sludge stream is held constant, which matches the usual baseline BSM1 operating story. |\n", + "| Reduced | Bioreactor layout | The original five biological compartments are lumped into one anoxic and one aerobic reactor. |\n", + "| Reduced | Biology and state space | The notebook keeps a reduced ASM1-inspired state vector and surrogate kinetics, so the reactor biology is illustrative rather than benchmark-identical. |\n", + "| Reduced | Clarifier | The ten-layer settler is replaced by an ideal solids separator so the example stays readable and fast. |\n", + "| Reduced | Influent scenario | The influent uses a synthetic disturbed dry-weather profile instead of the standard benchmark influent files. |\n", + "| Not included | Standard BSM1 KPIs | This notebook does not compute the official EQI or OCI indices, so the results are not directly benchmark-comparable. |" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "## Define the reduced BSM1 state and parameter set\n", + "\n", + "We start with a small immutable stream object and a parameter bundle. The stream tracks the reduced ASM1-style state variables that will move through the Plugboard process." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "STATE_FIELDS: tuple[str, ...] = (\n", + " \"s_s\",\n", + " \"x_s\",\n", + " \"x_bh\",\n", + " \"x_ba\",\n", + " \"s_o\",\n", + " \"s_no\",\n", + " \"s_nh\",\n", + " \"s_nd\",\n", + " \"x_nd\",\n", + ")\n", + "PARTICULATE_FIELDS: tuple[str, ...] = (\"x_s\", \"x_bh\", \"x_ba\", \"x_nd\")\n", + "EPS = 1e-6\n", + "\n", + "\n", + "class WastewaterStream(BaseModel):\n", + " \"\"\"Represents a simplified ASM1 stream used in the reduced BSM1 demo.\"\"\"\n", + "\n", + " model_config = ConfigDict(frozen=True)\n", + "\n", + " flow_m3_d: float\n", + " s_s: float\n", + " x_s: float\n", + " x_bh: float\n", + " x_ba: float\n", + " s_o: float\n", + " s_no: float\n", + " s_nh: float\n", + " s_nd: float\n", + " x_nd: float\n", + "\n", + " @property\n", + " def effluent_cod(self) -> float:\n", + " return self.s_s + self.x_s\n", + "\n", + " @property\n", + " def total_nitrogen(self) -> float:\n", + " return self.s_nh + self.s_no + self.s_nd + self.x_nd\n", + "\n", + "\n", + "class BSM1Parameters(BaseModel):\n", + " \"\"\"Parameter set for the reduced BSM1 demo.\"\"\"\n", + "\n", + " model_config = ConfigDict(frozen=True)\n", + "\n", + " dt_days: float = 5.0 / 1440.0\n", + " total_days: float = 14.0\n", + " v_anoxic_m3: float = 2000.0\n", + " v_aerobic_m3: float = 4000.0\n", + " q_in_avg_m3_d: float = 18446.0\n", + " q_return_m3_d: float = 18446.0\n", + " q_internal_bias_m3_d: float = 3.0 * 18446.0\n", + " q_internal_min_m3_d: float = 0.5 * 18446.0\n", + " q_internal_max_m3_d: float = 5.0 * 18446.0\n", + " oxygen_saturation_g_m3: float = 8.0\n", + " capture_efficiency: float = 0.995\n", + " do_setpoint_g_m3: float = 2.0\n", + " nitrate_setpoint_g_m3: float = 1.0\n", + " kla_bias_d_inv: float = 120.0\n", + " kla_min_d_inv: float = 30.0\n", + " kla_max_d_inv: float = 240.0\n", + " kla_rate_limit_d_inv_per_step: float = 20.0\n", + " do_kp: float = 27.5\n", + " do_ki: float = 110.0\n", + " nitrate_kp: float = 2500.0\n", + " nitrate_ki: float = 12000.0\n", + " y_h: float = 0.67\n", + " y_a: float = 0.24\n", + " i_xb: float = 0.08\n", + " mu_h_d_inv: float = 4.0\n", + " mu_a_d_inv: float = 0.5\n", + " b_h_d_inv: float = 0.3\n", + " b_a_d_inv: float = 0.05\n", + " k_s_g_m3: float = 10.0\n", + " k_oh_g_m3: float = 0.2\n", + " k_oa_g_m3: float = 0.4\n", + " k_no_g_m3: float = 0.5\n", + " k_nh_g_m3: float = 1.0\n", + " k_h_d_inv: float = 3.0\n", + " k_x_g_m3: float = 100.0\n", + " k_a_d_inv: float = 0.08\n", + " k_nd_g_m3: float = 1.0\n", + " eta_g: float = 0.8\n", + " eta_h: float = 0.4\n", + " decay_to_xs_fraction: float = 0.35\n", + " nitrification_oxygen_factor: float = 1.5\n", + "\n", + " @property\n", + " def total_steps(self) -> int:\n", + " return int(self.total_days / self.dt_days)" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## A plain-language guide to the wastewater variables\n", + "\n", + "If you do not work in wastewater treatment every day, the stream fields can look opaque. The reduced model keeps only the variables needed to tell the control story.\n", + "\n", + "- `flow_m3_d`: total stream flow in cubic metres per day.\n", + "- `s_s`: readily biodegradable soluble substrate. You can think of this as the easiest food for the heterotrophic biomass to consume.\n", + "- `x_s`: slowly biodegradable particulate substrate. This is still organic matter, but it must hydrolyse before the biomass can use it.\n", + "- `x_bh`: heterotrophic biomass, the organisms that mainly remove organic carbon and can denitrify when oxygen is scarce.\n", + "- `x_ba`: autotrophic biomass, the nitrifiers that oxidise ammonium when oxygen is available.\n", + "- `s_o`: dissolved oxygen, the main controlled variable in the aerobic reactor.\n", + "- `s_no`: oxidised nitrogen, represented here as nitrate and nitrite together. The internal recycle loop sends this back to the anoxic zone for denitrification.\n", + "- `s_nh`: ammonium nitrogen, the reduced nitrogen form that the aerobic biology tries to oxidise.\n", + "- `s_nd` and `x_nd`: soluble and particulate organic nitrogen pools. These provide a simple way to keep nitrogen tied to the organic matter balance.\n", + "\n", + "The notebook is still a reduced teaching model, but these variables are enough to connect the biology to the two feedback controllers and the main effluent trends." + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "## Helper functions for reactor balances and initial conditions\n", + "\n", + "The next cells keep the notebook readable by separating the low-level numerical helpers from the Plugboard components themselves." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "def _stream_kwargs(stream: WastewaterStream) -> dict[str, float]:\n", + " return {field: float(getattr(stream, field)) for field in STATE_FIELDS}\n", + "\n", + "\n", + "def _make_stream(flow_m3_d: float, values: dict[str, float]) -> WastewaterStream:\n", + " return WastewaterStream(flow_m3_d=flow_m3_d, **values)\n", + "\n", + "\n", + "def _non_negative_stream(stream: WastewaterStream) -> WastewaterStream:\n", + " return WastewaterStream(\n", + " flow_m3_d=max(stream.flow_m3_d, EPS),\n", + " **{field: max(getattr(stream, field), 0.0) for field in STATE_FIELDS},\n", + " )\n", + "\n", + "\n", + "def _with_flow(stream: WastewaterStream, flow_m3_d: float) -> WastewaterStream:\n", + " return WastewaterStream(flow_m3_d=max(flow_m3_d, EPS), **_stream_kwargs(stream))\n", + "\n", + "\n", + "def _blend_streams(streams: _t.Sequence[WastewaterStream]) -> WastewaterStream:\n", + " total_flow = sum(stream.flow_m3_d for stream in streams)\n", + " if total_flow <= EPS:\n", + " return WastewaterStream(\n", + " flow_m3_d=EPS,\n", + " s_s=0.0,\n", + " x_s=0.0,\n", + " x_bh=0.0,\n", + " x_ba=0.0,\n", + " s_o=0.0,\n", + " s_no=0.0,\n", + " s_nh=0.0,\n", + " s_nd=0.0,\n", + " x_nd=0.0,\n", + " )\n", + " blended = {\n", + " field: sum(stream.flow_m3_d * getattr(stream, field) for stream in streams) / total_flow\n", + " for field in STATE_FIELDS\n", + " }\n", + " return _make_stream(total_flow, blended)\n", + "\n", + "\n", + "def _advance_stream(\n", + " current: WastewaterStream,\n", + " inflows: _t.Sequence[WastewaterStream],\n", + " volume_m3: float,\n", + " params: BSM1Parameters,\n", + " kla_d_inv: float,\n", + ") -> WastewaterStream:\n", + " mixed = _blend_streams(inflows)\n", + " total_flow = mixed.flow_m3_d\n", + "\n", + " ss = max(current.s_s, 0.0)\n", + " xs = max(current.x_s, 0.0)\n", + " xbh = max(current.x_bh, 0.0)\n", + " xba = max(current.x_ba, 0.0)\n", + " so = max(current.s_o, 0.0)\n", + " sno = max(current.s_no, 0.0)\n", + " snh = max(current.s_nh, 0.0)\n", + " snd = max(current.s_nd, 0.0)\n", + " xnd = max(current.x_nd, 0.0)\n", + "\n", + " substrate_term = ss / (params.k_s_g_m3 + ss + EPS)\n", + " oxygen_term_h = so / (params.k_oh_g_m3 + so + EPS)\n", + " oxygen_term_a = so / (params.k_oa_g_m3 + so + EPS)\n", + " anoxic_term = params.k_oh_g_m3 / (params.k_oh_g_m3 + so + EPS)\n", + " nitrate_term = sno / (params.k_no_g_m3 + sno + EPS)\n", + " ammonium_term = snh / (params.k_nh_g_m3 + snh + EPS)\n", + " hydrolysis_term = xs / (params.k_x_g_m3 + xs + EPS)\n", + " hydrolysis_drive = oxygen_term_h + params.eta_h * anoxic_term * nitrate_term\n", + " organic_n_term = snd / (params.k_nd_g_m3 + snd + EPS)\n", + "\n", + " rho_aer_h = params.mu_h_d_inv * substrate_term * oxygen_term_h * xbh\n", + " rho_anox_h = (\n", + " params.eta_g * params.mu_h_d_inv * substrate_term * anoxic_term * nitrate_term * xbh\n", + " )\n", + " rho_auto = params.mu_a_d_inv * ammonium_term * oxygen_term_a * xba\n", + " rho_decay_h = params.b_h_d_inv * xbh\n", + " rho_decay_a = params.b_a_d_inv * xba\n", + " rho_ammonification = params.k_a_d_inv * organic_n_term * xbh\n", + " rho_hydrolysis = params.k_h_d_inv * hydrolysis_term * hydrolysis_drive * xbh\n", + " rho_hydrolysis_n = (xnd / max(xs, EPS)) * rho_hydrolysis if xs > EPS else 0.0\n", + "\n", + " reactions = {\n", + " \"s_s\": rho_hydrolysis - (rho_aer_h + rho_anox_h) / params.y_h,\n", + " \"x_s\": -rho_hydrolysis + params.decay_to_xs_fraction * (rho_decay_h + rho_decay_a),\n", + " \"x_bh\": rho_aer_h + rho_anox_h - rho_decay_h,\n", + " \"x_ba\": rho_auto - rho_decay_a,\n", + " \"s_o\": -((1.0 - params.y_h) / params.y_h) * rho_aer_h\n", + " - params.nitrification_oxygen_factor * rho_auto\n", + " + max(kla_d_inv, 0.0) * (params.oxygen_saturation_g_m3 - so),\n", + " \"s_no\": (1.0 / params.y_a) * rho_auto\n", + " - ((1.0 - params.y_h) / (2.86 * params.y_h)) * rho_anox_h,\n", + " \"s_nh\": -(params.i_xb * (rho_aer_h + rho_anox_h))\n", + " - ((1.0 / params.y_a) + params.i_xb) * rho_auto\n", + " + rho_ammonification\n", + " + params.i_xb * (rho_decay_h + rho_decay_a),\n", + " \"s_nd\": -rho_ammonification + rho_hydrolysis_n,\n", + " \"x_nd\": params.i_xb * (rho_decay_h + rho_decay_a) - rho_hydrolysis_n,\n", + " }\n", + "\n", + " dilution = total_flow / volume_m3\n", + " updated: dict[str, float] = {}\n", + " current_values = _stream_kwargs(current)\n", + " mixed_values = _stream_kwargs(mixed)\n", + " for field in STATE_FIELDS:\n", + " derivative = dilution * (mixed_values[field] - current_values[field]) + reactions[field]\n", + " updated[field] = current_values[field] + params.dt_days * derivative\n", + "\n", + " return _non_negative_stream(_make_stream(total_flow, updated))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "def _base_influent(params: BSM1Parameters) -> WastewaterStream:\n", + " return WastewaterStream(\n", + " flow_m3_d=params.q_in_avg_m3_d,\n", + " s_s=69.5,\n", + " x_s=202.32,\n", + " x_bh=28.17,\n", + " x_ba=0.0,\n", + " s_o=0.0,\n", + " s_no=0.0,\n", + " s_nh=31.56,\n", + " s_nd=6.95,\n", + " x_nd=10.59,\n", + " )\n", + "\n", + "\n", + "def _initial_anoxic_state(params: BSM1Parameters) -> WastewaterStream:\n", + " return WastewaterStream(\n", + " flow_m3_d=params.q_in_avg_m3_d + params.q_return_m3_d + params.q_internal_bias_m3_d,\n", + " s_s=28.0,\n", + " x_s=175.0,\n", + " x_bh=2250.0,\n", + " x_ba=110.0,\n", + " s_o=0.15,\n", + " s_no=1.1,\n", + " s_nh=16.0,\n", + " s_nd=3.4,\n", + " x_nd=11.0,\n", + " )\n", + "\n", + "\n", + "def _initial_aerobic_state(params: BSM1Parameters) -> WastewaterStream:\n", + " return WastewaterStream(\n", + " flow_m3_d=params.q_in_avg_m3_d + params.q_return_m3_d,\n", + " s_s=8.0,\n", + " x_s=110.0,\n", + " x_bh=2100.0,\n", + " x_ba=185.0,\n", + " s_o=params.do_setpoint_g_m3,\n", + " s_no=6.5,\n", + " s_nh=2.0,\n", + " s_nd=1.1,\n", + " x_nd=7.0,\n", + " )\n", + "\n", + "\n", + "def _initial_internal_recycle(params: BSM1Parameters) -> WastewaterStream:\n", + " return _with_flow(_initial_aerobic_state(params), params.q_internal_bias_m3_d)\n", + "\n", + "\n", + "def _initial_settler_streams(\n", + " params: BSM1Parameters,\n", + ") -> tuple[WastewaterStream, WastewaterStream, WastewaterStream]:\n", + " aerobic = _initial_aerobic_state(params)\n", + " feed = _with_flow(aerobic, params.q_in_avg_m3_d + params.q_return_m3_d)\n", + " eff_flow = max(feed.flow_m3_d - params.q_return_m3_d, EPS)\n", + " solids_to_return = params.capture_efficiency\n", + "\n", + " effluent_values = _stream_kwargs(feed)\n", + " return_values = _stream_kwargs(feed)\n", + " for field in PARTICULATE_FIELDS:\n", + " feed_mass = getattr(feed, field) * feed.flow_m3_d\n", + " effluent_values[field] = ((1.0 - solids_to_return) * feed_mass) / eff_flow\n", + " return_values[field] = (solids_to_return * feed_mass) / params.q_return_m3_d\n", + "\n", + " return (\n", + " feed,\n", + " _non_negative_stream(_make_stream(eff_flow, effluent_values)),\n", + " _non_negative_stream(_make_stream(params.q_return_m3_d, return_values)),\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "## Define the Plugboard components\n", + "\n", + "This is the key part of the example. Each major unit operation and control element becomes a normal Plugboard component, with explicit inputs and outputs that we will later wire together with connectors." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "class InfluentScenario(Component):\n", + " \"\"\"Generates a dynamic dry-weather-plus-disturbance influent profile.\"\"\"\n", + "\n", + " io = IO(\n", + " outputs=[\n", + " \"influent\",\n", + " \"time_day\",\n", + " \"do_setpoint\",\n", + " \"nitrate_setpoint\",\n", + " \"tick\",\n", + " \"influent_flow\",\n", + " \"influent_s_s\",\n", + " \"influent_s_nh\",\n", + " ]\n", + " )\n", + "\n", + " def __init__(self, params: BSM1Parameters, **kwargs: _t.Unpack[ComponentArgsDict]) -> None:\n", + " super().__init__(**kwargs)\n", + " self._params = params\n", + " self._step_index = 0\n", + " self._base = _base_influent(params)\n", + "\n", + " def _profile(self, time_day: float) -> WastewaterStream:\n", + " q_scale = 1.0 + 0.08 * math.sin(2.0 * math.pi * time_day)\n", + " q_scale += 0.04 * math.sin(4.0 * math.pi * time_day - 1.0)\n", + " cod_scale = 1.0 + 0.10 * math.sin(2.0 * math.pi * (time_day - 0.15))\n", + " nh_scale = 1.0 + 0.07 * math.sin(2.0 * math.pi * (time_day - 0.10))\n", + "\n", + " if 4.0 <= time_day < 6.0:\n", + " nh_scale += 0.25\n", + " if 8.0 <= time_day < 10.0:\n", + " q_scale += 0.35\n", + " cod_scale -= 0.10\n", + " nh_scale -= 0.05\n", + " if 11.0 <= time_day < 12.5:\n", + " cod_scale += 0.20\n", + " nh_scale += 0.10\n", + "\n", + " return WastewaterStream(\n", + " flow_m3_d=self._base.flow_m3_d * q_scale,\n", + " s_s=self._base.s_s * cod_scale,\n", + " x_s=self._base.x_s * cod_scale,\n", + " x_bh=self._base.x_bh,\n", + " x_ba=self._base.x_ba,\n", + " s_o=0.0,\n", + " s_no=0.0,\n", + " s_nh=self._base.s_nh * nh_scale,\n", + " s_nd=self._base.s_nd * nh_scale,\n", + " x_nd=self._base.x_nd * nh_scale,\n", + " )\n", + "\n", + " async def step(self) -> None:\n", + " if self._step_index >= self._params.total_steps:\n", + " await self.io.close()\n", + " return\n", + "\n", + " time_day = self._step_index * self._params.dt_days\n", + " self.influent = self._profile(time_day)\n", + " self.time_day = time_day\n", + " self.do_setpoint = self._params.do_setpoint_g_m3\n", + " self.nitrate_setpoint = self._params.nitrate_setpoint_g_m3\n", + " self.tick = self._step_index\n", + " self.influent_flow = self.influent.flow_m3_d\n", + " self.influent_s_s = self.influent.s_s\n", + " self.influent_s_nh = self.influent.s_nh\n", + "\n", + " self._step_index += 1\n", + "\n", + "\n", + "class PIController(Component):\n", + " \"\"\"Simple PI controller with output limits and anti-windup.\"\"\"\n", + "\n", + " io = IO(inputs=[\"measurement\", \"setpoint\", \"tick\"], outputs=[\"control_signal\"])\n", + "\n", + " def __init__(\n", + " self,\n", + " kp: float,\n", + " ki: float,\n", + " bias: float,\n", + " lower: float,\n", + " upper: float,\n", + " dt_days: float,\n", + " **kwargs: _t.Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._kp = kp\n", + " self._ki = ki\n", + " self._bias = bias\n", + " self._lower = lower\n", + " self._upper = upper\n", + " self._dt_days = dt_days\n", + " self._integral = 0.0\n", + "\n", + " async def step(self) -> None:\n", + " error = self.setpoint - self.measurement\n", + " proposed_integral = self._integral + self._ki * error * self._dt_days\n", + " raw_output = self._bias + (self._kp * error) + proposed_integral\n", + " bounded_output = min(max(raw_output, self._lower), self._upper)\n", + "\n", + " should_integrate = self._lower < raw_output < self._upper\n", + " should_integrate = should_integrate or (raw_output <= self._lower and error > 0.0)\n", + " should_integrate = should_integrate or (raw_output >= self._upper and error < 0.0)\n", + " if should_integrate:\n", + " self._integral = proposed_integral\n", + " raw_output = self._bias + (self._kp * error) + self._integral\n", + " bounded_output = min(max(raw_output, self._lower), self._upper)\n", + "\n", + " self.control_signal = bounded_output" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "class AnoxicReactor(Component):\n", + " \"\"\"Reduced anoxic bioreactor with ASM1-inspired nitrogen and COD dynamics.\"\"\"\n", + "\n", + " io = IO(\n", + " inputs=[\"influent\", \"internal_recycle\", \"return_sludge\", \"tick\"],\n", + " outputs=[\"mixed_liquor\", \"s_no\", \"anoxic_s_s\"],\n", + " )\n", + "\n", + " def __init__(\n", + " self,\n", + " params: BSM1Parameters,\n", + " initial_state: WastewaterStream,\n", + " **kwargs: _t.Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._params = params\n", + " self._stream_state = initial_state\n", + "\n", + " async def step(self) -> None:\n", + " self._stream_state = _advance_stream(\n", + " current=self._stream_state,\n", + " inflows=[self.influent, self.internal_recycle, self.return_sludge],\n", + " volume_m3=self._params.v_anoxic_m3,\n", + " params=self._params,\n", + " kla_d_inv=0.0,\n", + " )\n", + " self.mixed_liquor = self._stream_state\n", + " self.s_no = self._stream_state.s_no\n", + " self.anoxic_s_s = self._stream_state.s_s\n", + "\n", + "\n", + "class AerobicReactor(Component):\n", + " \"\"\"Reduced aerobic reactor with explicit aeration and internal recycle split.\"\"\"\n", + "\n", + " io = IO(\n", + " inputs=[\"feed\", \"recycle_flow\", \"kla\", \"tick\"],\n", + " outputs=[\n", + " \"to_settler\",\n", + " \"internal_recycle\",\n", + " \"s_o\",\n", + " \"aerobic_s_nh\",\n", + " \"internal_recycle_flow\",\n", + " \"applied_kla\",\n", + " ],\n", + " )\n", + "\n", + " def __init__(\n", + " self,\n", + " params: BSM1Parameters,\n", + " initial_state: WastewaterStream,\n", + " **kwargs: _t.Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._params = params\n", + " self._stream_state = initial_state\n", + " self._kla_state = params.kla_bias_d_inv\n", + "\n", + " async def step(self) -> None:\n", + " requested_kla = min(\n", + " max(self.kla, self._params.kla_min_d_inv),\n", + " self._params.kla_max_d_inv,\n", + " )\n", + " kla_delta = requested_kla - self._kla_state\n", + " max_delta = self._params.kla_rate_limit_d_inv_per_step\n", + " if kla_delta > max_delta:\n", + " self._kla_state += max_delta\n", + " elif kla_delta < -max_delta:\n", + " self._kla_state -= max_delta\n", + " else:\n", + " self._kla_state = requested_kla\n", + "\n", + " self._stream_state = _advance_stream(\n", + " current=self._stream_state,\n", + " inflows=[self.feed],\n", + " volume_m3=self._params.v_aerobic_m3,\n", + " params=self._params,\n", + " kla_d_inv=self._kla_state,\n", + " )\n", + "\n", + " recycle_flow = min(\n", + " max(self.recycle_flow, self._params.q_internal_min_m3_d),\n", + " self._params.q_internal_max_m3_d,\n", + " )\n", + " settler_flow = max(self.feed.flow_m3_d - recycle_flow, EPS)\n", + "\n", + " self.internal_recycle = _with_flow(self._stream_state, recycle_flow)\n", + " self.to_settler = _with_flow(self._stream_state, settler_flow)\n", + " self.s_o = self._stream_state.s_o\n", + " self.aerobic_s_nh = self._stream_state.s_nh\n", + " self.internal_recycle_flow = recycle_flow\n", + " self.applied_kla = self._kla_state\n", + "\n", + "\n", + "class IdealSettler(Component):\n", + " \"\"\"Ideal solids separator used as a pedagogical replacement for the 10-layer clarifier.\"\"\"\n", + "\n", + " io = IO(\n", + " inputs=[\"feed\", \"tick\"],\n", + " outputs=[\n", + " \"return_sludge\",\n", + " \"effluent_cod\",\n", + " \"effluent_s_s\",\n", + " \"effluent_s_nh\",\n", + " \"effluent_s_no\",\n", + " \"effluent_total_n\",\n", + " \"return_sludge_x_bh\",\n", + " ],\n", + " )\n", + "\n", + " def __init__(self, params: BSM1Parameters, **kwargs: _t.Unpack[ComponentArgsDict]) -> None:\n", + " super().__init__(**kwargs)\n", + " self._params = params\n", + "\n", + " async def step(self) -> None:\n", + " effluent_flow = max(self.feed.flow_m3_d - self._params.q_return_m3_d, EPS)\n", + " effluent_values = _stream_kwargs(self.feed)\n", + " return_values = _stream_kwargs(self.feed)\n", + "\n", + " for field in PARTICULATE_FIELDS:\n", + " feed_mass = getattr(self.feed, field) * self.feed.flow_m3_d\n", + " effluent_values[field] = (\n", + " (1.0 - self._params.capture_efficiency) * feed_mass / effluent_flow\n", + " )\n", + " return_values[field] = (\n", + " self._params.capture_efficiency * feed_mass / self._params.q_return_m3_d\n", + " )\n", + "\n", + " effluent = _non_negative_stream(_make_stream(effluent_flow, effluent_values))\n", + " self.return_sludge = _non_negative_stream(\n", + " _make_stream(self._params.q_return_m3_d, return_values)\n", + " )\n", + " self.effluent_cod = effluent.effluent_cod\n", + " self.effluent_s_s = effluent.s_s\n", + " self.effluent_s_nh = effluent.s_nh\n", + " self.effluent_s_no = effluent.s_no\n", + " self.effluent_total_n = effluent.total_nitrogen\n", + " self.return_sludge_x_bh = self.return_sludge.x_bh" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "## Assemble the Plugboard process\n", + "\n", + "\n", + "\n", + "Once the components are defined, the rest is standard Plugboard work: instantiate the components, create connectors, and build a LocalProcess. The writer is connected directly to scalar outputs from the process components, so a few extra outputs exist only to make file writing and plotting explicit inside the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "def build_process(\n", + " output_path: str | Path = \"bsm1_results.csv\",\n", + " params: BSM1Parameters | None = None,\n", + ") -> LocalProcess:\n", + " \"\"\"Builds the reduced BSM1 demo process.\"\"\"\n", + "\n", + " params = params or BSM1Parameters()\n", + " initial_anoxic = _initial_anoxic_state(params)\n", + " initial_aerobic = _initial_aerobic_state(params)\n", + " initial_internal = _initial_internal_recycle(params)\n", + " initial_settler_feed, _, initial_return = _initial_settler_streams(params)\n", + "\n", + " def connect(source: str, target: str) -> AsyncioConnector:\n", + " return AsyncioConnector(spec=ConnectorSpec(source=source, target=target))\n", + "\n", + " influent = InfluentScenario(name=\"influent\", params=params)\n", + " nitrate_controller = PIController(\n", + " name=\"nitrate_controller\",\n", + " kp=params.nitrate_kp,\n", + " ki=params.nitrate_ki,\n", + " bias=params.q_internal_bias_m3_d,\n", + " lower=params.q_internal_min_m3_d,\n", + " upper=params.q_internal_max_m3_d,\n", + " dt_days=params.dt_days,\n", + " initial_values={\n", + " \"measurement\": [initial_anoxic.s_no],\n", + " \"setpoint\": [params.nitrate_setpoint_g_m3],\n", + " },\n", + " )\n", + " do_controller = PIController(\n", + " name=\"do_controller\",\n", + " kp=params.do_kp,\n", + " ki=params.do_ki,\n", + " bias=params.kla_bias_d_inv,\n", + " lower=params.kla_min_d_inv,\n", + " upper=params.kla_max_d_inv,\n", + " dt_days=params.dt_days,\n", + " initial_values={\n", + " \"measurement\": [initial_aerobic.s_o],\n", + " \"setpoint\": [params.do_setpoint_g_m3],\n", + " },\n", + " )\n", + " anoxic = AnoxicReactor(\n", + " name=\"anoxic_reactor\",\n", + " params=params,\n", + " initial_state=initial_anoxic,\n", + " initial_values={\n", + " \"internal_recycle\": [initial_internal],\n", + " \"return_sludge\": [initial_return],\n", + " },\n", + " )\n", + " aerobic = AerobicReactor(\n", + " name=\"aerobic_reactor\",\n", + " params=params,\n", + " initial_state=initial_aerobic,\n", + " initial_values={\n", + " \"feed\": [initial_anoxic],\n", + " \"recycle_flow\": [params.q_internal_bias_m3_d],\n", + " \"kla\": [params.kla_bias_d_inv],\n", + " },\n", + " )\n", + " settler = IdealSettler(\n", + " name=\"ideal_settler\",\n", + " params=params,\n", + " initial_values={\"feed\": [initial_settler_feed]},\n", + " )\n", + " writer = FileWriter(\n", + " name=\"writer\",\n", + " path=str(output_path),\n", + " field_names=[\n", + " \"time_day_out\",\n", + " \"influent_flow\",\n", + " \"influent_s_s\",\n", + " \"influent_s_nh\",\n", + " \"anoxic_s_no\",\n", + " \"anoxic_s_s\",\n", + " \"aerobic_s_o\",\n", + " \"aerobic_s_nh\",\n", + " \"effluent_cod\",\n", + " \"effluent_s_s\",\n", + " \"effluent_s_nh\",\n", + " \"effluent_s_no\",\n", + " \"effluent_total_n\",\n", + " \"return_sludge_x_bh\",\n", + " \"internal_recycle_flow\",\n", + " \"kla_out\",\n", + " \"requested_kla_out\",\n", + " \"do_setpoint_out\",\n", + " \"nitrate_setpoint_out\",\n", + " ],\n", + " )\n", + "\n", + " components = [\n", + " influent,\n", + " nitrate_controller,\n", + " do_controller,\n", + " anoxic,\n", + " aerobic,\n", + " settler,\n", + " writer,\n", + " ]\n", + " connectors = [\n", + " connect(\"influent.influent\", \"anoxic_reactor.influent\"),\n", + " connect(\"influent.tick\", \"nitrate_controller.tick\"),\n", + " connect(\"influent.tick\", \"do_controller.tick\"),\n", + " connect(\"influent.tick\", \"anoxic_reactor.tick\"),\n", + " connect(\"influent.tick\", \"aerobic_reactor.tick\"),\n", + " connect(\"influent.tick\", \"ideal_settler.tick\"),\n", + " connect(\"influent.time_day\", \"writer.time_day_out\"),\n", + " connect(\"influent.influent_flow\", \"writer.influent_flow\"),\n", + " connect(\"influent.influent_s_s\", \"writer.influent_s_s\"),\n", + " connect(\"influent.influent_s_nh\", \"writer.influent_s_nh\"),\n", + " connect(\"influent.do_setpoint\", \"do_controller.setpoint\"),\n", + " connect(\"influent.do_setpoint\", \"writer.do_setpoint_out\"),\n", + " connect(\"influent.nitrate_setpoint\", \"nitrate_controller.setpoint\"),\n", + " connect(\"influent.nitrate_setpoint\", \"writer.nitrate_setpoint_out\"),\n", + " connect(\"anoxic_reactor.mixed_liquor\", \"aerobic_reactor.feed\"),\n", + " connect(\"anoxic_reactor.s_no\", \"nitrate_controller.measurement\"),\n", + " connect(\"anoxic_reactor.s_no\", \"writer.anoxic_s_no\"),\n", + " connect(\"anoxic_reactor.anoxic_s_s\", \"writer.anoxic_s_s\"),\n", + " connect(\"nitrate_controller.control_signal\", \"aerobic_reactor.recycle_flow\"),\n", + " connect(\"aerobic_reactor.internal_recycle\", \"anoxic_reactor.internal_recycle\"),\n", + " connect(\"aerobic_reactor.to_settler\", \"ideal_settler.feed\"),\n", + " connect(\"aerobic_reactor.s_o\", \"do_controller.measurement\"),\n", + " connect(\"aerobic_reactor.s_o\", \"writer.aerobic_s_o\"),\n", + " connect(\"aerobic_reactor.aerobic_s_nh\", \"writer.aerobic_s_nh\"),\n", + " connect(\"aerobic_reactor.internal_recycle_flow\", \"writer.internal_recycle_flow\"),\n", + " connect(\"aerobic_reactor.applied_kla\", \"writer.kla_out\"),\n", + " connect(\"do_controller.control_signal\", \"aerobic_reactor.kla\"),\n", + " connect(\"do_controller.control_signal\", \"writer.requested_kla_out\"),\n", + " connect(\"ideal_settler.return_sludge\", \"anoxic_reactor.return_sludge\"),\n", + " connect(\"ideal_settler.effluent_cod\", \"writer.effluent_cod\"),\n", + " connect(\"ideal_settler.effluent_s_s\", \"writer.effluent_s_s\"),\n", + " connect(\"ideal_settler.effluent_s_nh\", \"writer.effluent_s_nh\"),\n", + " connect(\"ideal_settler.effluent_s_no\", \"writer.effluent_s_no\"),\n", + " connect(\"ideal_settler.effluent_total_n\", \"writer.effluent_total_n\"),\n", + " connect(\"ideal_settler.return_sludge_x_bh\", \"writer.return_sludge_x_bh\"),\n", + " ]\n", + "\n", + " return LocalProcess(components=components, connectors=connectors)\n", + "\n", + "\n", + "params = BSM1Parameters()\n", + "process = build_process(output_path=\"bsm1_results.csv\", params=params)" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "## Diagram of the process\n", + "\n", + "\n", + "\n", + "The process graph below is generated from the actual LocalProcess wiring. That makes it a useful visual check that the internal recycle, return sludge, controllers, and direct writer connections are connected in the way the code intends.\n", + "\n", + "\n", + "\n", + "![Process Diagram](https://mermaid.ink/img/pako:eNrtmcFqwzAMhl_F-BjaFwijbAwGgR1GethlEFxbpQbHGopLByXvPlq3Y2EsjCGzkPn-_z8fsiwZfJQaDchSbh0e9E5REI_1ixfC-q3bgw-3R9Ht1CuUgnDvDZiFcGoDrhTVRbHW4BVZvNnQqiiuvqIQvVguV0J5fLO6IVA6IH0fd3fW1VEWs4bWUyIrGYFGMjDC9ESooevqizDmXG2z5_E2kArQaPSB0LlRsur-QxXjvrrZAQ3-km1gZMeabscrINz8CC0Kh2xDMzucNaBc00EIo6dZnWTrqLokfTbO_lZOtOnz8PqbgrEOm2RDgp9ysuVL1XgpctOdN3Nmmr2ahJKrMVOw8a7Yf1c-3tXLl5bsGjMjcr0R2Lo40XOdl2-aVWOiYog54xzIhrGUB-vg-ayJAVGfKTJFpsgUmSJTZIpMMV8KuZAtUKuskeVRhh20p08fA1u1d0H2_TsDnDKC)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "diagram_url = MermaidDiagram.from_process(process).url\n", + "print(diagram_url)\n", + "print(markdown_diagram(process))" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "## Inspect the parameter set\n", + "\n", + "The values below are the main operating and controller parameters exposed by the demo model. The dissolved-oxygen controller gains are intentionally softer than the earlier notebook revision, and the aerobic reactor now applies a small slew limit to `KLa` so the aeration actuator cannot jump unrealistically from one 5-minute sample to the next. That keeps the DO trace easier to read while preserving the basic control story." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "parameter_table = pd.DataFrame(\n", + " [\n", + " (\"Anoxic volume\", params.v_anoxic_m3, \"m^3\"),\n", + " (\"Aerobic volume\", params.v_aerobic_m3, \"m^3\"),\n", + " (\"Average influent flow\", params.q_in_avg_m3_d, \"m^3/d\"),\n", + " (\"Return sludge flow\", params.q_return_m3_d, \"m^3/d\"),\n", + " (\"Internal recycle bias\", params.q_internal_bias_m3_d, \"m^3/d\"),\n", + " (\"DO setpoint\", params.do_setpoint_g_m3, \"g/m^3\"),\n", + " (\"Nitrate setpoint\", params.nitrate_setpoint_g_m3, \"g N/m^3\"),\n", + " (\"KLa bias\", params.kla_bias_d_inv, \"d^-1\"),\n", + " (\"KLa slew limit per step\", params.kla_rate_limit_d_inv_per_step, \"d^-1\"),\n", + " (\"DO controller Kp\", params.do_kp, \"-\"),\n", + " (\"DO controller Ki\", params.do_ki, \"d^-1\"),\n", + " (\"Nitrate controller Kp\", params.nitrate_kp, \"-\"),\n", + " (\"Nitrate controller Ki\", params.nitrate_ki, \"d^-1\"),\n", + " (\"Time step\", params.dt_days * 24 * 60, \"minutes\"),\n", + " (\"Simulation horizon\", params.total_days, \"days\"),\n", + " ],\n", + " columns=[\"Parameter\", \"Value\", \"Unit\"],\n", + ")\n", + "parameter_table" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "## Run the simulation\n", + "\n", + "For the run cell we create a fresh process instance so the execution state is clean even if you re-run the notebook sections out of order." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "results_path = Path(\"bsm1_results.csv\")\n", + "results_path.unlink(missing_ok=True)\n", + "process = build_process(output_path=results_path, params=params)\n", + "async with process:\n", + " await process.run()\n", + "df = pd.read_csv(results_path)\n", + "df.head()" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "## Summarize the closed-loop response\n", + "\n", + "\n", + "\n", + "Official BSM1 studies often report standardized quality and cost indices. This demo keeps the summary lighter: it looks at average effluent quality and average actuator usage over the final three days of the run.\n", + "\n", + "\n", + "\n", + "These numbers are tutorial-scale surrogates rather than benchmark-comparable BSM1 KPIs. They are most useful for checking whether the simplified control story behaves plausibly under the imposed disturbances." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23", + "metadata": {}, + "outputs": [], + "source": [ + "final_window = df[df[\"time_day_out\"] >= (params.total_days - 3.0)]\n", + "summary = pd.DataFrame(\n", + " {\n", + " \"Metric\": [\n", + " \"Average effluent COD surrogate (S_S + X_S)\",\n", + " \"Average effluent soluble substrate\",\n", + " \"Average effluent ammonium\",\n", + " \"Average effluent nitrate\",\n", + " \"Average effluent total nitrogen\",\n", + " \"Average internal recycle flow\",\n", + " \"Average KLa\",\n", + " ],\n", + " \"Value\": [\n", + " final_window[\"effluent_cod\"].mean(),\n", + " final_window[\"effluent_s_s\"].mean(),\n", + " final_window[\"effluent_s_nh\"].mean(),\n", + " final_window[\"effluent_s_no\"].mean(),\n", + " final_window[\"effluent_total_n\"].mean(),\n", + " final_window[\"internal_recycle_flow\"].mean(),\n", + " final_window[\"kla_out\"].mean(),\n", + " ],\n", + " }\n", + ")\n", + "summary" + ] + }, + { + "cell_type": "markdown", + "id": "24", + "metadata": {}, + "source": [ + "## Visualize plant behaviour\n", + "\n", + "The first panel shows the imposed influent disturbances. The middle panels show the controlled states and the controller moves. The last panel tracks the effluent quality variables that are most useful for reading this reduced model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25", + "metadata": {}, + "outputs": [], + "source": [ + "fig = make_subplots(\n", + " rows=4,\n", + " cols=1,\n", + " shared_xaxes=True,\n", + " vertical_spacing=0.05,\n", + " subplot_titles=(\n", + " \"Influent disturbances\",\n", + " \"Controlled reactor states\",\n", + " \"Manipulated variables\",\n", + " \"Effluent quality\",\n", + " ),\n", + ")\n", + "\n", + "fig.add_trace(\n", + " go.Scatter(x=df[\"time_day_out\"], y=df[\"influent_flow\"], name=\"Influent flow\"),\n", + " row=1,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(x=df[\"time_day_out\"], y=df[\"influent_s_nh\"], name=\"Influent ammonium\"),\n", + " row=1,\n", + " col=1,\n", + ")\n", + "\n", + "fig.add_trace(\n", + " go.Scatter(x=df[\"time_day_out\"], y=df[\"aerobic_s_o\"], name=\"Aerobic DO\"),\n", + " row=2,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df[\"time_day_out\"],\n", + " y=df[\"do_setpoint_out\"],\n", + " name=\"DO setpoint\",\n", + " line={\"dash\": \"dash\"},\n", + " ),\n", + " row=2,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(x=df[\"time_day_out\"], y=df[\"anoxic_s_no\"], name=\"Anoxic nitrate\"),\n", + " row=2,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df[\"time_day_out\"],\n", + " y=df[\"nitrate_setpoint_out\"],\n", + " name=\"Nitrate setpoint\",\n", + " line={\"dash\": \"dot\"},\n", + " ),\n", + " row=2,\n", + " col=1,\n", + ")\n", + "\n", + "fig.add_trace(\n", + " go.Scatter(x=df[\"time_day_out\"], y=df[\"requested_kla_out\"], name=\"Requested KLa\"),\n", + " row=3,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df[\"time_day_out\"],\n", + " y=df[\"kla_out\"],\n", + " name=\"Applied KLa\",\n", + " line={\"width\": 3},\n", + " ),\n", + " row=3,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df[\"time_day_out\"],\n", + " y=df[\"internal_recycle_flow\"],\n", + " name=\"Internal recycle flow\",\n", + " ),\n", + " row=3,\n", + " col=1,\n", + ")\n", + "\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df[\"time_day_out\"],\n", + " y=df[\"effluent_cod\"],\n", + " name=\"Effluent COD surrogate\",\n", + " ),\n", + " row=4,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df[\"time_day_out\"],\n", + " y=df[\"effluent_s_nh\"],\n", + " name=\"Effluent ammonium\",\n", + " ),\n", + " row=4,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df[\"time_day_out\"],\n", + " y=df[\"effluent_total_n\"],\n", + " name=\"Effluent total N\",\n", + " ),\n", + " row=4,\n", + " col=1,\n", + ")\n", + "\n", + "fig.update_layout(height=1100, width=950, template=\"plotly_white\", legend={\"orientation\": \"h\"})\n", + "fig.update_xaxes(title_text=\"Time [days]\", row=4, col=1)\n", + "fig.update_yaxes(title_text=\"Influent\", row=1, col=1)\n", + "fig.update_yaxes(title_text=\"Controlled states\", row=2, col=1)\n", + "fig.update_yaxes(title_text=\"Controller outputs\", row=3, col=1)\n", + "fig.update_yaxes(title_text=\"Effluent quality\", row=4, col=1)\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "26", + "metadata": {}, + "source": [ + "## Interpretation\n", + "\n", + "A few things are worth checking when you run the notebook:\n", + "\n", + "- the internal recycle should move as anoxic nitrate drifts away from its setpoint\n", + "- the requested `KLa` can still move quickly, but the applied `KLa` should ramp more smoothly because the aerobic reactor now enforces a per-step slew limit\n", + "- the aerobic DO trace should still react to disturbances, but with less sample-to-sample chatter than the earlier direct-actuator version\n", + "- effluent ammonium and the COD surrogate should worsen during the imposed load increases, then recover as the controllers respond\n", + "- the ideal settler keeps solids handling simple, so this notebook emphasizes biological control and process wiring rather than clarifier hydraulics\n", + "\n", + "This version combines the earlier DO gain reduction with a simple actuator model. The controller still decides how much aeration it wants, but the plant applies that request through a limited-rate `KLa` change, which is closer to how a real blower or valve train behaves than an instantaneous jump. That makes the plotted trajectories easier to read without changing the tutorial scope of the model.\n", + "\n", + "Treat these trends as qualitative checks on the reduced control topology, not as benchmark-faithful BSM1 performance claims. Natural extensions are adding a more detailed clarifier, splitting the bioreactor back into multiple compartments, or computing official BSM1-style quality and operating-cost indices." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "plugboard (3.12.8)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/demos/transport/.meta.yml b/examples/demos/transport/.meta.yml new file mode 100644 index 00000000..fed0f44a --- /dev/null +++ b/examples/demos/transport/.meta.yml @@ -0,0 +1,2 @@ +tags: + - transport diff --git a/examples/demos/transport/001_traffic_simulation/.meta.yml b/examples/demos/transport/001_traffic_simulation/.meta.yml new file mode 100644 index 00000000..cd359cb7 --- /dev/null +++ b/examples/demos/transport/001_traffic_simulation/.meta.yml @@ -0,0 +1,2 @@ +tags: + - events diff --git a/examples/demos/transport/001_traffic_simulation/traffic-simulation.ipynb b/examples/demos/transport/001_traffic_simulation/traffic-simulation.ipynb new file mode 100644 index 00000000..770a9bba --- /dev/null +++ b/examples/demos/transport/001_traffic_simulation/traffic-simulation.ipynb @@ -0,0 +1,1045 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Traffic intersection simulation\n", + "\n", + "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/plugboard-dev/plugboard/blob/main/examples/demos/transport/001_traffic_simulation/traffic-simulation.ipynb)\n", + "\n", + "[![Open in GitHub Codespaces](https://github.com/codespaces/badge.svg)](https://codespaces.new/plugboard-dev/plugboard)\n", + "\n", + "***Event-driven simulations are non-deterministic and not suitable for running in distributed environments, because Plugboard won't guarantee the order in which events get processed.***\n", + "\n", + "This demo simulates a dual-intersection traffic corridor using Plugboard's **event-driven** architecture.\n", + "Vehicles arrive at **Junction A**, wait for a green signal, cross the intersection, travel along a\n", + "connecting road with stochastic travel time, then repeat the process at **Junction B** before exiting the system.\n", + "\n", + "### Components\n", + "\n", + "| Component | Role |\n", + "|---|---|\n", + "| **SimulationClock** | Emits integer time steps (1 step = 1 second) |\n", + "| **VehicleSource** | Generates vehicles via a Poisson process |\n", + "| **TrafficLight** ×2 | FSM controllers cycling GREEN → YELLOW → RED |\n", + "| **Junction** ×2 | FIFO queues that release vehicles when the signal is green |\n", + "| **ConnectingRoad** | Adds triangular-distributed travel time between junctions |\n", + "| **MetricsCollector** | Subscribes to *all* event types for post-simulation analysis |\n", + "\n", + "### Event types\n", + "\n", + "| Event | Trigger |\n", + "|---|---|\n", + "| `VehicleArrivedEvent` | A vehicle enters a junction queue |\n", + "| `VehicleDepartedEvent` | A vehicle passes through a junction |\n", + "| `SignalChangedEvent` | A traffic light transitions between states |\n", + "| `CongestionAlertEvent` | A junction's queue exceeds its congestion threshold |\n", + "\n", + "Regular connectors carry continuous data (clock time, signal state, queue length), while\n", + "**events** capture discrete occurrences — the two paradigms work together in a single process.\n", + "\n", + "We first run a **synchronised** baseline (both lights share the same phase), then show how\n", + "offsetting Junction B's green phase to create a **green wave** improves throughput." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "# Install plugboard and dependencies for Google Colab\n", + "!pip install -q plugboard plotly" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import random\n", + "\n", + "import pandas as pd\n", + "import plotly.graph_objects as go\n", + "from plotly.subplots import make_subplots\n", + "\n", + "from plugboard.connector import AsyncioConnector\n", + "from plugboard.process import LocalProcess\n", + "from plugboard.schemas import ConnectorSpec\n", + "\n", + "connect = lambda src, tgt: AsyncioConnector(spec=ConnectorSpec(source=src, target=tgt))" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "## Component definitions\n", + "\n", + "The simulation is built from six Plugboard components. Each one declares its I/O via the class-level\n", + "`io` attribute, which Plugboard uses to wire up both regular field connectors and event routes.\n", + "\n", + "Four lightweight event types carry discrete signals around the system:\n", + "\n", + "| Event | Payload |\n", + "|---|---|\n", + "| `VehicleArrivedEvent` | `vehicle_id`, `target_junction`, `time` |\n", + "| `VehicleDepartedEvent` | `vehicle_id`, `from_junction`, `wait_time`, `time` |\n", + "| `SignalChangedEvent` | `junction`, `new_state`, `old_state`, `time` |\n", + "| `CongestionAlertEvent` | `junction`, `queue_length`, `threshold`, `time` |" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "from collections import deque\n", + "from typing import Any, ClassVar, Unpack\n", + "\n", + "from pydantic import BaseModel\n", + "\n", + "from plugboard.component import Component, IOController as IO\n", + "from plugboard.events import Event\n", + "from plugboard.schemas import ComponentArgsDict\n", + "\n", + "\n", + "# ── Event data models ────────────────────────────────────────────────────────\n", + "\n", + "\n", + "class VehicleArrivalData(BaseModel):\n", + " \"\"\"Payload for a vehicle arrival event.\"\"\"\n", + "\n", + " vehicle_id: int\n", + " target_junction: str\n", + " time: float\n", + "\n", + "\n", + "class VehicleDepartureData(BaseModel):\n", + " \"\"\"Payload for a vehicle departure event.\"\"\"\n", + "\n", + " vehicle_id: int\n", + " from_junction: str\n", + " wait_time: float\n", + " time: float\n", + "\n", + "\n", + "class SignalChangeData(BaseModel):\n", + " \"\"\"Payload for a traffic signal change event.\"\"\"\n", + "\n", + " junction: str\n", + " new_state: str\n", + " old_state: str\n", + " time: float\n", + "\n", + "\n", + "class CongestionData(BaseModel):\n", + " \"\"\"Payload for a congestion alert event.\"\"\"\n", + "\n", + " junction: str\n", + " queue_length: int\n", + " threshold: int\n", + " time: float\n", + "\n", + "\n", + "# ── Events ───────────────────────────────────────────────────────────────────\n", + "\n", + "\n", + "class VehicleArrivedEvent(Event):\n", + " \"\"\"Emitted when a vehicle arrives at a junction queue.\"\"\"\n", + "\n", + " type: ClassVar[str] = \"vehicle_arrived\"\n", + " data: VehicleArrivalData\n", + "\n", + "\n", + "class VehicleDepartedEvent(Event):\n", + " \"\"\"Emitted when a vehicle passes through a junction.\"\"\"\n", + "\n", + " type: ClassVar[str] = \"vehicle_departed\"\n", + " data: VehicleDepartureData\n", + "\n", + "\n", + "class SignalChangedEvent(Event):\n", + " \"\"\"Emitted when a traffic light changes state.\"\"\"\n", + "\n", + " type: ClassVar[str] = \"signal_changed\"\n", + " data: SignalChangeData\n", + "\n", + "\n", + "class CongestionAlertEvent(Event):\n", + " \"\"\"Emitted when a junction queue exceeds the congestion threshold.\"\"\"\n", + "\n", + " type: ClassVar[str] = \"congestion_alert\"\n", + " data: CongestionData\n", + "\n", + "\n", + "# ── Components ───────────────────────────────────────────────────────────────\n", + "\n", + "\n", + "class SimulationClock(Component):\n", + " \"\"\"Drives the simulation by emitting integer time steps (1 step = 1 second).\"\"\"\n", + "\n", + " io = IO(outputs=[\"time\"])\n", + "\n", + " def __init__(\n", + " self,\n", + " duration: int = 600,\n", + " **kwargs: Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._duration: int = duration\n", + " self._time: int = 0\n", + "\n", + " async def step(self) -> None:\n", + " self.time = self._time\n", + " self._time += 1\n", + " if self._time > self._duration:\n", + " self._logger.info(\"Simulation complete\", duration=self._duration)\n", + " await self.io.close()\n", + "\n", + "\n", + "class VehicleSource(Component):\n", + " \"\"\"Generates vehicles with Poisson arrivals.\n", + "\n", + " At each step a vehicle arrives with probability ``arrival_rate``.\n", + " Each vehicle is assigned a unique ID and a ``VehicleArrivedEvent``\n", + " is emitted targeting the specified junction.\n", + " \"\"\"\n", + "\n", + " io = IO(inputs=[\"time\"], output_events=[VehicleArrivedEvent])\n", + "\n", + " def __init__(\n", + " self,\n", + " arrival_rate: float = 0.4,\n", + " target_junction: str = \"A\",\n", + " **kwargs: Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._arrival_rate: float = arrival_rate\n", + " self._target_junction: str = target_junction\n", + " self._next_id: int = 1\n", + "\n", + " async def step(self) -> None:\n", + " if random.random() < self._arrival_rate: # noqa: S311\n", + " vid = self._next_id\n", + " self._next_id += 1\n", + " self.io.queue_event(\n", + " VehicleArrivedEvent(\n", + " source=self.name,\n", + " data=VehicleArrivalData(\n", + " vehicle_id=vid,\n", + " target_junction=self._target_junction,\n", + " time=float(self.time),\n", + " ),\n", + " )\n", + " )\n", + "\n", + "\n", + "class TrafficLight(Component):\n", + " \"\"\"FSM-based traffic signal controller.\n", + "\n", + " Cycles through GREEN → YELLOW → RED with configurable phase durations.\n", + " The ``offset`` parameter shifts the starting phase so that two lights\n", + " can be synchronised or staggered to create green-wave effects.\n", + " A ``SignalChangedEvent`` is emitted on every state transition.\n", + " \"\"\"\n", + "\n", + " io = IO(\n", + " inputs=[\"time\"],\n", + " outputs=[\"signal\"],\n", + " output_events=[SignalChangedEvent],\n", + " )\n", + "\n", + " def __init__(\n", + " self,\n", + " junction: str = \"A\",\n", + " green_duration: int = 30,\n", + " yellow_duration: int = 5,\n", + " red_duration: int = 30,\n", + " offset: int = 0,\n", + " **kwargs: Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._junction: str = junction\n", + " self._green: int = green_duration\n", + " self._yellow: int = yellow_duration\n", + " self._red: int = red_duration\n", + " self._offset: int = offset\n", + " self._cycle: int = green_duration + yellow_duration + red_duration\n", + " self._signal_state: str | None = None\n", + "\n", + " def _compute_state(self, t: int) -> str:\n", + " phase = (t + self._offset) % self._cycle\n", + " if phase < self._green:\n", + " return \"green\"\n", + " if phase < self._green + self._yellow:\n", + " return \"yellow\"\n", + " return \"red\"\n", + "\n", + " async def step(self) -> None:\n", + " t = int(self.time)\n", + " new = self._compute_state(t)\n", + " if new != self._signal_state:\n", + " old = self._signal_state or \"init\"\n", + " self.io.queue_event(\n", + " SignalChangedEvent(\n", + " source=self.name,\n", + " data=SignalChangeData(\n", + " junction=self._junction,\n", + " new_state=new,\n", + " old_state=old,\n", + " time=float(t),\n", + " ),\n", + " )\n", + " )\n", + " self._signal_state = new\n", + " self.signal = self._signal_state\n", + "\n", + "\n", + "class Junction(Component):\n", + " \"\"\"FIFO vehicle queue at a signalised intersection.\n", + "\n", + " Vehicles arrive via ``VehicleArrivedEvent`` and are filtered by\n", + " ``junction_name``. When the signal is green, vehicles are released\n", + " at ``service_rate`` (probability per step). Each release produces a\n", + " ``VehicleDepartedEvent``. If the queue exceeds\n", + " ``congestion_threshold`` a ``CongestionAlertEvent`` is emitted.\n", + " \"\"\"\n", + "\n", + " io = IO(\n", + " inputs=[\"time\", \"signal\"],\n", + " outputs=[\"queue_length\", \"total_passed\"],\n", + " input_events=[VehicleArrivedEvent],\n", + " output_events=[VehicleDepartedEvent, CongestionAlertEvent],\n", + " )\n", + "\n", + " def __init__(\n", + " self,\n", + " junction_name: str = \"A\",\n", + " service_rate: float = 0.8,\n", + " capacity: int = 100,\n", + " congestion_threshold: int = 20,\n", + " **kwargs: Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._junction_name: str = junction_name\n", + " self._service_rate: float = service_rate\n", + " self._capacity: int = capacity\n", + " self._congestion_threshold: int = congestion_threshold\n", + " self._queue: deque[tuple[int, float]] = deque()\n", + " self._total_passed: int = 0\n", + " self._congestion_active: bool = False\n", + "\n", + " @VehicleArrivedEvent.handler\n", + " async def on_vehicle_arrived(self, event: VehicleArrivedEvent) -> None:\n", + " \"\"\"Enqueue vehicles targeted at this junction.\"\"\"\n", + " if event.data.target_junction == self._junction_name:\n", + " if len(self._queue) < self._capacity:\n", + " self._queue.append((event.data.vehicle_id, event.data.time))\n", + "\n", + " async def step(self) -> None:\n", + " t = float(self.time)\n", + " if self.signal == \"green\" and self._queue:\n", + " if random.random() < self._service_rate: # noqa: S311\n", + " vehicle_id, arrival_time = self._queue.popleft()\n", + " self._total_passed += 1\n", + " self.io.queue_event(\n", + " VehicleDepartedEvent(\n", + " source=self.name,\n", + " data=VehicleDepartureData(\n", + " vehicle_id=vehicle_id,\n", + " from_junction=self._junction_name,\n", + " wait_time=t - arrival_time,\n", + " time=t,\n", + " ),\n", + " )\n", + " )\n", + " qlen = len(self._queue)\n", + " if qlen > self._congestion_threshold and not self._congestion_active:\n", + " self._congestion_active = True\n", + " self._logger.warning(\n", + " \"Congestion detected\",\n", + " junction=self._junction_name,\n", + " queue_length=qlen,\n", + " )\n", + " self.io.queue_event(\n", + " CongestionAlertEvent(\n", + " source=self.name,\n", + " data=CongestionData(\n", + " junction=self._junction_name,\n", + " queue_length=qlen,\n", + " threshold=self._congestion_threshold,\n", + " time=t,\n", + " ),\n", + " )\n", + " )\n", + " elif qlen <= self._congestion_threshold:\n", + " self._congestion_active = False\n", + " self.queue_length = qlen\n", + " self.total_passed = self._total_passed\n", + "\n", + "\n", + "class ConnectingRoad(Component):\n", + " \"\"\"Road segment between two junctions with stochastic travel time.\n", + "\n", + " Handles ``VehicleDepartedEvent`` from the upstream junction, buffers\n", + " the vehicle for a random travel time drawn from a triangular\n", + " distribution, then emits a ``VehicleArrivedEvent`` at the downstream\n", + " junction.\n", + " \"\"\"\n", + "\n", + " io = IO(\n", + " inputs=[\"time\"],\n", + " input_events=[VehicleDepartedEvent],\n", + " output_events=[VehicleArrivedEvent],\n", + " )\n", + "\n", + " def __init__(\n", + " self,\n", + " from_junction: str = \"A\",\n", + " to_junction: str = \"B\",\n", + " min_travel: float = 10.0,\n", + " mode_travel: float = 15.0,\n", + " max_travel: float = 25.0,\n", + " **kwargs: Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._from: str = from_junction\n", + " self._to: str = to_junction\n", + " self._min: float = min_travel\n", + " self._mode: float = mode_travel\n", + " self._max: float = max_travel\n", + " self._in_transit: list[tuple[int, float]] = []\n", + "\n", + " @VehicleDepartedEvent.handler\n", + " async def on_departed(self, event: VehicleDepartedEvent) -> None:\n", + " \"\"\"Buffer a departing vehicle with a random travel delay.\"\"\"\n", + " if event.data.from_junction == self._from:\n", + " travel = random.triangular(self._min, self._max, self._mode) # noqa: S311\n", + " self._in_transit.append((event.data.vehicle_id, event.data.time + travel))\n", + "\n", + " async def step(self) -> None:\n", + " t = float(self.time)\n", + " remaining: list[tuple[int, float]] = []\n", + " for vid, release in self._in_transit:\n", + " if t >= release:\n", + " self.io.queue_event(\n", + " VehicleArrivedEvent(\n", + " source=self.name,\n", + " data=VehicleArrivalData(\n", + " vehicle_id=vid,\n", + " target_junction=self._to,\n", + " time=t,\n", + " ),\n", + " )\n", + " )\n", + " else:\n", + " remaining.append((vid, release))\n", + " self._in_transit = remaining\n", + "\n", + "\n", + "class MetricsCollector(Component):\n", + " \"\"\"Subscribes to all simulation events and records them for analysis.\n", + "\n", + " After the simulation finishes, access the ``*_log`` lists and\n", + " ``queue_history`` to build plots and compute KPIs.\n", + " \"\"\"\n", + "\n", + " io = IO(\n", + " inputs=[\"time\", \"queue_a\", \"queue_b\", \"passed_a\", \"passed_b\"],\n", + " input_events=[\n", + " VehicleArrivedEvent,\n", + " VehicleDepartedEvent,\n", + " CongestionAlertEvent,\n", + " SignalChangedEvent,\n", + " ],\n", + " )\n", + "\n", + " def __init__(self, **kwargs: Unpack[ComponentArgsDict]) -> None:\n", + " super().__init__(**kwargs)\n", + " self.arrival_log: list[dict[str, Any]] = []\n", + " self.departure_log: list[dict[str, Any]] = []\n", + " self.congestion_log: list[dict[str, Any]] = []\n", + " self.signal_log: list[dict[str, Any]] = []\n", + " self.queue_history: list[dict[str, Any]] = []\n", + "\n", + " @VehicleArrivedEvent.handler\n", + " async def on_arrival(self, event: VehicleArrivedEvent) -> None:\n", + " \"\"\"Record vehicle arrival.\"\"\"\n", + " self.arrival_log.append(event.data.model_dump())\n", + "\n", + " @VehicleDepartedEvent.handler\n", + " async def on_departure(self, event: VehicleDepartedEvent) -> None:\n", + " \"\"\"Record vehicle departure.\"\"\"\n", + " self.departure_log.append(event.data.model_dump())\n", + "\n", + " @CongestionAlertEvent.handler\n", + " async def on_congestion(self, event: CongestionAlertEvent) -> None:\n", + " \"\"\"Record congestion alert.\"\"\"\n", + " self.congestion_log.append(event.data.model_dump())\n", + "\n", + " @SignalChangedEvent.handler\n", + " async def on_signal(self, event: SignalChangedEvent) -> None:\n", + " \"\"\"Record signal change.\"\"\"\n", + " self.signal_log.append(event.data.model_dump())\n", + "\n", + " async def step(self) -> None:\n", + " self.queue_history.append(\n", + " {\n", + " \"time\": float(self.time),\n", + " \"queue_a\": int(self.queue_a),\n", + " \"queue_b\": int(self.queue_b),\n", + " \"passed_a\": int(self.passed_a),\n", + " \"passed_b\": int(self.passed_b),\n", + " }\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "## Assembling the process\n", + "\n", + "We create 8 components and wire them together:\n", + "\n", + "- **Field connectors** carry the clock time to every component, signal states from lights to junctions,\n", + " and queue metrics into the collector.\n", + "- **Event connectors** are built automatically from the `input_events` / `output_events` declared on each\n", + " component's `IOController`. Plugboard routes every emitted event to all components that subscribe to its type.\n", + "\n", + "The baseline scenario uses synchronised signals — both junctions share the same GREEN/YELLOW/RED phase." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "random.seed(42) # reproducible results\n", + "\n", + "DURATION = 600 # seconds (10 minutes)\n", + "ARRIVAL_RATE = 0.4 # vehicles per second (~24 vehicles/min)\n", + "GREEN = 30 # seconds\n", + "YELLOW = 5\n", + "RED = 30\n", + "\n", + "components = [\n", + " SimulationClock(name=\"clock\", duration=DURATION),\n", + " VehicleSource(name=\"source\", arrival_rate=ARRIVAL_RATE, target_junction=\"A\"),\n", + " TrafficLight(\n", + " name=\"light-a\", junction=\"A\", green_duration=GREEN, yellow_duration=YELLOW, red_duration=RED\n", + " ),\n", + " TrafficLight(\n", + " name=\"light-b\", junction=\"B\", green_duration=GREEN, yellow_duration=YELLOW, red_duration=RED\n", + " ),\n", + " Junction(name=\"junction-a\", junction_name=\"A\", congestion_threshold=15),\n", + " Junction(name=\"junction-b\", junction_name=\"B\", congestion_threshold=15),\n", + " ConnectingRoad(name=\"road-ab\", from_junction=\"A\", to_junction=\"B\"),\n", + " MetricsCollector(name=\"metrics\"),\n", + "]\n", + "\n", + "connectors = [\n", + " # Clock drives every component\n", + " connect(\"clock.time\", \"source.time\"),\n", + " connect(\"clock.time\", \"light-a.time\"),\n", + " connect(\"clock.time\", \"light-b.time\"),\n", + " connect(\"clock.time\", \"junction-a.time\"),\n", + " connect(\"clock.time\", \"junction-b.time\"),\n", + " connect(\"clock.time\", \"road-ab.time\"),\n", + " connect(\"clock.time\", \"metrics.time\"),\n", + " # Signal state into junctions\n", + " connect(\"light-a.signal\", \"junction-a.signal\"),\n", + " connect(\"light-b.signal\", \"junction-b.signal\"),\n", + " # Queue metrics into collector\n", + " connect(\"junction-a.queue_length\", \"metrics.queue_a\"),\n", + " connect(\"junction-b.queue_length\", \"metrics.queue_b\"),\n", + " connect(\"junction-a.total_passed\", \"metrics.passed_a\"),\n", + " connect(\"junction-b.total_passed\", \"metrics.passed_b\"),\n", + "]\n", + "\n", + "# Event connectors are built automatically from IO declarations\n", + "event_connectors = AsyncioConnector.builder().build_event_connectors(components)\n", + "\n", + "process = LocalProcess(components=components, connectors=connectors + event_connectors)\n", + "print(f\"{len(process.components)} components, {len(process.connectors)} connectors\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "from plugboard.diagram import MermaidDiagram\n", + "\n", + "diagram_url = MermaidDiagram.from_process(process).url\n", + "print(diagram_url)" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "In the process diagram dotted lines represent the event-driven portions of the model.\n", + "\n", + "![Process Diagram](https://mermaid.ink/img/pako:eNrtmF1v2yAUhv-KxWWUL3_EX5qqVWlvpvammXYxRYqwwTYbhgjjrF2U_z4ISbtOS6cGcre7QN7z8vhw8LG9BSVHGOSgovxH2UAhvbuHJfO8kvLy-8et1zVwjXNP8J4hjIYehQWmubcgbU-hJJzNtfBDIa4Gg33MYODtvNHoyut4L0p82uILbkhJ8WIvMwYmRDu4QaCkbuQInvb4LGBVkfJO64zBIcQ1Q_F-hsIdw7eelVrwVio-HTQm_iXiAhTFuykc5kJwiEbwDYQ5ZwyrZVn9oKTG4hDkjqLFUpCyO-1xbwRzTqmC4cKYHMKOHHb17aw07ErcWW3YXorDjbG9nEug_M_Kbygb03xWUAiywejFrcGPf3aoa6O53WAmNcTY0bFxxmBxZo4MCK9V7_8HxM1B9IrCyf2Usxp3mmwFKRby7xTzZ9W1Fr2icFMVHakZpCv1HKRWOpGLxV4zN5ILMNhWhVuK8-vC0Y7YPEQakPMT6qLFaoJzy8pFb7Vd3759aAKbA-6Owqao3fQvp7mworDNhe1t3_ZsgiFosWghQerldauJlkA2uMVLkKufCFewp3IJlkxLYS_54omVIJeix0OgYOsG5BWknRr1awQlviGwFrA9SszkLSLq1nSco4ofq9EWyKe1fmuuSSeVvdrSitR6vhdUTTdSrrt8MtF_j2sim74Yl7yddATpV-xmk8WTOIhTGIQ4TkI4C0NUFn6WVkHkVyiZ-gEEu90QrCH7yvkzkxrqRR5B7vvxOE3iKJn5WTabRkE6BE8gT7JxlKZJEiSzMAijKFIeP_cG03Hip2E2TeIgnc6yMFUBeH9x9-YDwP47wBDUQmf0kCWsdlTM1cZKkM92vwAU_pKm)" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "## Run the baseline simulation\n", + "\n", + "With both signals synchronised, vehicles released from Junction A arrive at Junction B\n", + "during its RED phase — they queue up and wait." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "async with process:\n", + " await process.run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "mc = process.components[\"metrics\"]\n", + "\n", + "df_queue = pd.DataFrame(mc.queue_history)\n", + "df_arrivals = pd.DataFrame(mc.arrival_log)\n", + "df_departures = pd.DataFrame(mc.departure_log)\n", + "df_signals = pd.DataFrame(mc.signal_log)\n", + "df_congestion = pd.DataFrame(mc.congestion_log)\n", + "\n", + "total_arrived = len(df_arrivals[df_arrivals[\"target_junction\"] == \"A\"])\n", + "passed_a = df_queue[\"passed_a\"].iloc[-1]\n", + "passed_b = df_queue[\"passed_b\"].iloc[-1]\n", + "\n", + "print(f\"Vehicles generated : {total_arrived}\")\n", + "print(f\"Through Junction A : {passed_a}\")\n", + "print(f\"Through Junction B : {passed_b} (exited system)\")\n", + "print(f\"System throughput : {passed_b / total_arrived:.1%}\")\n", + "print(f\"Congestion alerts : {len(df_congestion)}\")\n", + "if len(df_departures):\n", + " print(\n", + " f\"Mean wait time (A) : {df_departures[df_departures['from_junction'] == 'A']['wait_time'].mean():.1f}s\"\n", + " )\n", + " print(\n", + " f\"Mean wait time (B) : {df_departures[df_departures['from_junction'] == 'B']['wait_time'].mean():.1f}s\"\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "### Queue dynamics\n", + "\n", + "The chart below shows how queue lengths at each junction oscillate with the signal cycle.\n", + "Green and red background bands indicate the signal phase at Junction A." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "CYCLE = GREEN + YELLOW + RED\n", + "\n", + "fig = make_subplots(\n", + " rows=2,\n", + " cols=1,\n", + " shared_xaxes=True,\n", + " vertical_spacing=0.08,\n", + " subplot_titles=(\"Junction A queue\", \"Junction B queue\"),\n", + ")\n", + "\n", + "# Signal-state bands for Junction A\n", + "for start in range(0, DURATION, CYCLE):\n", + " fig.add_vrect(\n", + " x0=start, x1=start + GREEN, fillcolor=\"green\", opacity=0.08, line_width=0, row=1, col=1\n", + " )\n", + " fig.add_vrect(\n", + " x0=start + GREEN + YELLOW,\n", + " x1=start + CYCLE,\n", + " fillcolor=\"red\",\n", + " opacity=0.08,\n", + " line_width=0,\n", + " row=1,\n", + " col=1,\n", + " )\n", + "\n", + "# Signal-state bands for Junction B (no offset in baseline)\n", + "for start in range(0, DURATION, CYCLE):\n", + " fig.add_vrect(\n", + " x0=start, x1=start + GREEN, fillcolor=\"green\", opacity=0.08, line_width=0, row=2, col=1\n", + " )\n", + " fig.add_vrect(\n", + " x0=start + GREEN + YELLOW,\n", + " x1=start + CYCLE,\n", + " fillcolor=\"red\",\n", + " opacity=0.08,\n", + " line_width=0,\n", + " row=2,\n", + " col=1,\n", + " )\n", + "\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_queue[\"time\"], y=df_queue[\"queue_a\"], name=\"Queue A\", line=dict(color=\"royalblue\")\n", + " ),\n", + " row=1,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(x=df_queue[\"time\"], y=df_queue[\"queue_b\"], name=\"Queue B\", line=dict(color=\"coral\")),\n", + " row=2,\n", + " col=1,\n", + ")\n", + "\n", + "# Mark congestion alerts\n", + "if len(df_congestion):\n", + " for _, row in df_congestion.iterrows():\n", + " r = 1 if row[\"junction\"] == \"A\" else 2\n", + " fig.add_vline(x=row[\"time\"], line_dash=\"dot\", line_color=\"red\", opacity=0.6, row=r, col=1)\n", + "\n", + "fig.update_layout(\n", + " height=500,\n", + " template=\"plotly_white\",\n", + " showlegend=False,\n", + " title_text=\"Baseline: synchronised signals\",\n", + ")\n", + "fig.update_xaxes(title_text=\"Time (s)\", row=2, col=1)\n", + "fig.update_yaxes(title_text=\"Vehicles\", row=1, col=1)\n", + "fig.update_yaxes(title_text=\"Vehicles\", row=2, col=1)\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "### Wait-time distributions\n", + "\n", + "Vehicles arriving during a red phase wait longer. The distribution at Junction B is typically\n", + "wider because it includes both signal delays and the connecting-road variance." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "fig = go.Figure()\n", + "for junc, colour in [(\"A\", \"royalblue\"), (\"B\", \"coral\")]:\n", + " waits = df_departures[df_departures[\"from_junction\"] == junc][\"wait_time\"]\n", + " if len(waits):\n", + " fig.add_trace(\n", + " go.Histogram(\n", + " x=waits, name=f\"Junction {junc}\", opacity=0.7, marker_color=colour, nbinsx=30\n", + " )\n", + " )\n", + "\n", + "fig.update_layout(\n", + " barmode=\"overlay\",\n", + " template=\"plotly_white\",\n", + " title=\"Vehicle wait-time distribution\",\n", + " xaxis_title=\"Wait time (s)\",\n", + " yaxis_title=\"Count\",\n", + ")\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "## Experiment: green-wave offset\n", + "\n", + "When both signals share the same phase, vehicles released from Junction A during its green\n", + "phase arrive at Junction B roughly 10–25 seconds later — often hitting B's red phase.\n", + "\n", + "A **green wave** staggers Junction B's phase so that its green window aligns with the arrival\n", + "of vehicles from A. The optimal offset is approximately:\n", + "\n", + "$$\\text{offset}_B = \\text{cycle length} - \\text{avg. travel time} \\approx 65 - 15 = 50\\text{s}$$\n", + "\n", + "Let's re-run the simulation with this offset and compare." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "def build_process(light_b_offset: int = 0, seed: int = 42) -> LocalProcess:\n", + " \"\"\"Create a fresh traffic-corridor process with a given offset for Junction B.\"\"\"\n", + " random.seed(seed)\n", + " comps = [\n", + " SimulationClock(name=\"clock\", duration=DURATION),\n", + " VehicleSource(name=\"source\", arrival_rate=ARRIVAL_RATE, target_junction=\"A\"),\n", + " TrafficLight(\n", + " name=\"light-a\",\n", + " junction=\"A\",\n", + " green_duration=GREEN,\n", + " yellow_duration=YELLOW,\n", + " red_duration=RED,\n", + " ),\n", + " TrafficLight(\n", + " name=\"light-b\",\n", + " junction=\"B\",\n", + " green_duration=GREEN,\n", + " yellow_duration=YELLOW,\n", + " red_duration=RED,\n", + " offset=light_b_offset,\n", + " ),\n", + " Junction(name=\"junction-a\", junction_name=\"A\", congestion_threshold=15),\n", + " Junction(name=\"junction-b\", junction_name=\"B\", congestion_threshold=15),\n", + " ConnectingRoad(name=\"road-ab\", from_junction=\"A\", to_junction=\"B\"),\n", + " MetricsCollector(name=\"metrics\"),\n", + " ]\n", + " conns = [\n", + " connect(\"clock.time\", \"source.time\"),\n", + " connect(\"clock.time\", \"light-a.time\"),\n", + " connect(\"clock.time\", \"light-b.time\"),\n", + " connect(\"clock.time\", \"junction-a.time\"),\n", + " connect(\"clock.time\", \"junction-b.time\"),\n", + " connect(\"clock.time\", \"road-ab.time\"),\n", + " connect(\"clock.time\", \"metrics.time\"),\n", + " connect(\"light-a.signal\", \"junction-a.signal\"),\n", + " connect(\"light-b.signal\", \"junction-b.signal\"),\n", + " connect(\"junction-a.queue_length\", \"metrics.queue_a\"),\n", + " connect(\"junction-b.queue_length\", \"metrics.queue_b\"),\n", + " connect(\"junction-a.total_passed\", \"metrics.passed_a\"),\n", + " connect(\"junction-b.total_passed\", \"metrics.passed_b\"),\n", + " ]\n", + " evt = AsyncioConnector.builder().build_event_connectors(comps)\n", + " return LocalProcess(components=comps, connectors=conns + evt)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "OFFSET = CYCLE - 15 # ≈ mode travel time\n", + "\n", + "proc_wave = build_process(light_b_offset=OFFSET)\n", + "async with proc_wave:\n", + " await proc_wave.run()\n", + "\n", + "mc_wave = proc_wave.components[\"metrics\"]\n", + "df_queue_wave = pd.DataFrame(mc_wave.queue_history)\n", + "df_dep_wave = pd.DataFrame(mc_wave.departure_log)\n", + "df_cong_wave = pd.DataFrame(mc_wave.congestion_log)\n", + "\n", + "print(f\"Green-wave offset: {OFFSET}s\")\n", + "print(f\"Through Junction B: {df_queue_wave['passed_b'].iloc[-1]} (was {passed_b})\")\n", + "print(f\"Congestion alerts : {len(df_cong_wave)} (was {len(df_congestion)})\")\n", + "if len(df_dep_wave):\n", + " print(\n", + " f\"Mean wait (B) : {df_dep_wave[df_dep_wave['from_junction'] == 'B']['wait_time'].mean():.1f}s\"\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "### Side-by-side comparison\n", + "\n", + "The top row shows the **synchronised** baseline; the bottom row shows the **green-wave** scenario.\n", + "Notice how Junction B's queue is dramatically reduced when its green phase is aligned with\n", + "arriving traffic from Junction A." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "fig = make_subplots(\n", + " rows=2,\n", + " cols=2,\n", + " shared_xaxes=True,\n", + " vertical_spacing=0.12,\n", + " horizontal_spacing=0.08,\n", + " subplot_titles=(\n", + " \"Baseline — Junction A\",\n", + " \"Baseline — Junction B\",\n", + " \"Green wave — Junction A\",\n", + " \"Green wave — Junction B\",\n", + " ),\n", + ")\n", + "\n", + "# Baseline\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_queue[\"time\"], y=df_queue[\"queue_a\"], line=dict(color=\"royalblue\"), showlegend=False\n", + " ),\n", + " row=1,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_queue[\"time\"], y=df_queue[\"queue_b\"], line=dict(color=\"coral\"), showlegend=False\n", + " ),\n", + " row=1,\n", + " col=2,\n", + ")\n", + "\n", + "# Green wave\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_queue_wave[\"time\"],\n", + " y=df_queue_wave[\"queue_a\"],\n", + " line=dict(color=\"royalblue\"),\n", + " showlegend=False,\n", + " ),\n", + " row=2,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_queue_wave[\"time\"],\n", + " y=df_queue_wave[\"queue_b\"],\n", + " line=dict(color=\"coral\"),\n", + " showlegend=False,\n", + " ),\n", + " row=2,\n", + " col=2,\n", + ")\n", + "\n", + "for r in [1, 2]:\n", + " for c in [1, 2]:\n", + " fig.update_yaxes(title_text=\"Queue\", row=r, col=c)\n", + "fig.update_xaxes(title_text=\"Time (s)\", row=2, col=1)\n", + "fig.update_xaxes(title_text=\"Time (s)\", row=2, col=2)\n", + "\n", + "fig.update_layout(\n", + " height=500,\n", + " template=\"plotly_white\",\n", + " title_text=\"Queue length comparison: synchronised vs green wave\",\n", + ")\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "fig = go.Figure()\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_queue[\"time\"],\n", + " y=df_queue[\"passed_b\"],\n", + " name=\"Baseline\",\n", + " line=dict(color=\"coral\", dash=\"dash\"),\n", + " )\n", + ")\n", + "fig.add_trace(\n", + " go.Scatter(\n", + " x=df_queue_wave[\"time\"],\n", + " y=df_queue_wave[\"passed_b\"],\n", + " name=\"Green wave\",\n", + " line=dict(color=\"seagreen\"),\n", + " )\n", + ")\n", + "fig.update_layout(\n", + " template=\"plotly_white\",\n", + " title=\"Cumulative vehicles through Junction B\",\n", + " xaxis_title=\"Time (s)\",\n", + " yaxis_title=\"Vehicles exited\",\n", + ")\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "This demo illustrated several strengths of Plugboard's event-driven architecture applied to transport simulation:\n", + "\n", + "1. **Hybrid communication** — regular connectors carry continuous state (time, signal phase, queue lengths)\n", + " while discrete events (arrivals, departures, congestion alerts) propagate through the system.\n", + "2. **FSM logic** — traffic lights use a finite state machine that emits `SignalChangedEvent` only on transitions,\n", + " avoiding unnecessary work when the state is unchanged.\n", + "3. **Stochastic processes** — vehicle arrivals follow a Poisson process and road travel times use a triangular\n", + " distribution, giving the simulation realistic variability.\n", + "4. **Event monitoring** — the `MetricsCollector` subscribes to every event type, enabling rich post-simulation\n", + " analysis without coupling the analytics to the simulation logic.\n", + "5. **What-if analysis** — by wrapping process creation in a helper function, we can easily sweep parameters\n", + " (like signal offset) and compare scenarios.\n", + "\n", + "Try experimenting with different arrival rates, signal timings, or road travel times to explore\n", + "how the system responds under stress." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "plugboard (3.12.8)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/demos/transport/002_bus_terminal/.meta.yml b/examples/demos/transport/002_bus_terminal/.meta.yml new file mode 100644 index 00000000..cd359cb7 --- /dev/null +++ b/examples/demos/transport/002_bus_terminal/.meta.yml @@ -0,0 +1,2 @@ +tags: + - events diff --git a/examples/demos/transport/002_bus_terminal/bus-terminal.ipynb b/examples/demos/transport/002_bus_terminal/bus-terminal.ipynb new file mode 100644 index 00000000..30f6bff5 --- /dev/null +++ b/examples/demos/transport/002_bus_terminal/bus-terminal.ipynb @@ -0,0 +1,1179 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Bus Bunching Simulation\n", + "\n", + "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/plugboard-dev/plugboard/blob/main/examples/demos/transport/001-bus-terminal/bus-terminal.ipynb)\n", + "\n", + "[![Open in GitHub Codespaces](https://github.com/codespaces/badge.svg)](https://codespaces.new/plugboard-dev/plugboard)\n", + "\n", + "***Event-driven simulations are non-deterministic and not suitable for running in distributed environments, because Plugboard won't guarantee the order in which events get processed.***\n", + "\n", + "This notebook demonstrates how Plugboard's **event-driven simulation** naturally models the feedback dynamics that cause **bus bunching** — one of the most pervasive problems in public transit.\n", + "\n", + "**Bus bunching** occurs when a delayed bus picks up more waiting passengers, increasing its dwell time and delaying it further, while the following bus finds fewer passengers and speeds ahead. This positive feedback loop causes buses to cluster together, ruining service regularity.\n", + "\n", + "We model a transit corridor with staggered bus departures:\n", + "- **Bus stops** that accumulate passengers via a Poisson process\n", + "- **Buses** that travel between stops, dwell to board/alight passengers, and experience traffic signal delays\n", + "- **A central dispatch controller** that implements a headway-based holding strategy to mitigate bunching\n", + "\n", + "We run the simulation twice — once **without control** to observe bunching, then **with holding control** — and compare the results." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## Install and import dependencies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "# Install plugboard for Google Colab\n", + "%pip install -q plugboard plotly" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "import typing as _t\n", + "\n", + "import numpy as np\n", + "import plotly.express as px\n", + "import plotly.graph_objects as go\n", + "from plotly.subplots import make_subplots\n", + "from pydantic import BaseModel\n", + "\n", + "from plugboard.component import Component, IOController as IO\n", + "from plugboard.connector import AsyncioConnector\n", + "from plugboard.diagram import MermaidDiagram\n", + "from plugboard.events import Event\n", + "from plugboard.process import LocalProcess\n", + "from plugboard.schemas import ComponentArgsDict, ConnectorSpec\n", + "\n", + "connect = lambda src, tgt: AsyncioConnector(spec=ConnectorSpec(source=src, target=tgt))" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## Route configuration\n", + "\n", + "We model a single bus corridor with 15 stops. Buses depart from stop 0 at regular headways and travel one-way to the final stop.\n", + "\n", + "Key parameters:\n", + "- **15 stops** with 400–600m spacing (~7.5 km total corridor)\n", + "- **Bus speed**: 8 m/s (~29 km/h, typical urban bus)\n", + "- **6 buses** dispatched at even headways from the first stop\n", + "- **Capacity**: 80 passengers per bus\n", + "- **Passenger arrival rates** vary by stop — higher in the central corridor\n", + "- **Simulation timestep**: 10 seconds of real time per tick\n", + "- **Initial perturbation**: bus-0 suffers an extra 120s delay at its 3rd stop (e.g. wheelchair boarding), which seeds the bunching cascade" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "# --- Route configuration ---\n", + "NUM_STOPS = 15\n", + "NUM_BUSES = 6\n", + "BUS_CAPACITY = 80\n", + "BUS_SPEED = 8.0 # m/s (~29 km/h)\n", + "TICK_SECONDS = 10 # each simulation tick = 10 real seconds\n", + "\n", + "# Distance between consecutive stops (meters)\n", + "STOP_DISTANCES = [500, 450, 400, 550, 500, 600, 450, 500, 550, 400, 500, 450, 500, 550]\n", + "\n", + "# Cumulative stop positions along the route\n", + "STOP_POSITIONS = [0.0]\n", + "for d in STOP_DISTANCES:\n", + " STOP_POSITIONS.append(STOP_POSITIONS[-1] + d)\n", + "ROUTE_LENGTH = STOP_POSITIONS[-1]\n", + "\n", + "# Headway between bus departures from stop 0 (in ticks)\n", + "DEPARTURE_HEADWAY_TICKS = 18 # 180 seconds = 3 minutes\n", + "\n", + "# Total simulation ticks: enough for all buses to complete the route\n", + "TRAVEL_TIME_TICKS = math.ceil(ROUTE_LENGTH / (BUS_SPEED * TICK_SECONDS))\n", + "TOTAL_TICKS = (\n", + " (NUM_BUSES - 1) * DEPARTURE_HEADWAY_TICKS + TRAVEL_TIME_TICKS + 250\n", + ") # extra margin for dwell + holding\n", + "\n", + "# Base passenger arrival rate per stop per tick (Poisson lambda)\n", + "# Higher at central stops, creating more dwell time variation\n", + "BASE_ARRIVAL_RATES = [0.4, 0.6, 1.0, 1.4, 1.8, 2.2, 2.5, 2.5, 2.2, 1.8, 1.4, 1.0, 0.6, 0.4, 0.2]\n", + "\n", + "# Alighting probability at each stop\n", + "ALIGHT_PROBS = [\n", + " 0.02,\n", + " 0.03,\n", + " 0.05,\n", + " 0.06,\n", + " 0.08,\n", + " 0.10,\n", + " 0.12,\n", + " 0.12,\n", + " 0.10,\n", + " 0.08,\n", + " 0.08,\n", + " 0.10,\n", + " 0.12,\n", + " 0.15,\n", + " 0.30,\n", + "]\n", + "\n", + "# Dwell time parameters\n", + "DWELL_FIXED = 8.0 # fixed dead time (door open/close), seconds\n", + "DWELL_PER_PAX = 2.5 # seconds per boarding/alighting passenger\n", + "\n", + "# Traffic signal parameters\n", + "NUM_SIGNALS = 7\n", + "SIGNAL_GREEN_FRAC = 0.50\n", + "\n", + "# Perturbation: extra delay on bus-0 at stop 2 (simulates an incident)\n", + "PERTURBATION_BUS = \"bus-0\"\n", + "PERTURBATION_STOP = 2\n", + "PERTURBATION_TICKS = 12 # 120 seconds extra delay\n", + "\n", + "# Headway control threshold (seconds)\n", + "HEADWAY_THRESHOLD = 30.0\n", + "\n", + "print(f\"Route: {ROUTE_LENGTH:.0f}m, {NUM_STOPS} stops\")\n", + "print(f\"Fleet: {NUM_BUSES} buses at {DEPARTURE_HEADWAY_TICKS * TICK_SECONDS}s headway\")\n", + "print(f\"Simulation: {TOTAL_TICKS} ticks = {TOTAL_TICKS * TICK_SECONDS / 60:.0f} min\")" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "## Define events\n", + "\n", + "Events model the discrete interactions between buses, stops, and the dispatch controller:\n", + "\n", + "| Event | Direction | Purpose |\n", + "|-------|-----------|---------|\n", + "| `BusArrivedAtStop` | Bus → Stop | Triggers passenger boarding/alighting |\n", + "| `PassengersBoarded` | Stop → Bus | Reports boarding results and dwell time |\n", + "| `BusPositionUpdate` | Bus → Dispatch | Periodic position report for headway tracking |\n", + "| `HoldBus` | Dispatch → Bus | Instruction to hold at a stop |" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "# --- Event data models ---\n", + "\n", + "\n", + "class BusArrivedData(BaseModel):\n", + " bus_id: str\n", + " stop_index: int\n", + " current_load: int\n", + " tick: int\n", + "\n", + "\n", + "class PassengersBoardedData(BaseModel):\n", + " bus_id: str\n", + " stop_index: int\n", + " boarded: int\n", + " alighted: int\n", + " new_load: int\n", + " dwell_ticks: int\n", + "\n", + "\n", + "class BusPositionData(BaseModel):\n", + " bus_id: str\n", + " position: float\n", + " stop_index: int\n", + " load: int\n", + " tick: int\n", + " phase: str\n", + "\n", + "\n", + "class HoldBusData(BaseModel):\n", + " bus_id: str\n", + " hold_ticks: int\n", + "\n", + "\n", + "# --- Event classes ---\n", + "\n", + "\n", + "class BusArrivedAtStop(Event):\n", + " type: _t.ClassVar[str] = \"bus_arrived_at_stop\"\n", + " data: BusArrivedData\n", + "\n", + "\n", + "class PassengersBoarded(Event):\n", + " type: _t.ClassVar[str] = \"passengers_boarded\"\n", + " data: PassengersBoardedData\n", + "\n", + "\n", + "class BusPositionUpdate(Event):\n", + " type: _t.ClassVar[str] = \"bus_position_update\"\n", + " data: BusPositionData\n", + "\n", + "\n", + "class HoldBus(Event):\n", + " type: _t.ClassVar[str] = \"hold_bus\"\n", + " data: HoldBusData" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "## Create the Bus Stop component\n", + "\n", + "Each `BusStop` maintains a passenger queue. On every tick, new passengers arrive via a Poisson process. When a bus arrives (via `BusArrivedAtStop`), the stop:\n", + "\n", + "1. Computes alighting passengers (fraction of onboard load)\n", + "2. Computes boarding passengers (limited by bus capacity)\n", + "3. Calculates dwell time: $T_d = t_0 + t_{pax} \\cdot (n_{board} + n_{alight})$\n", + "4. Records the arrival for headway analysis\n", + "5. Returns a `PassengersBoarded` event\n", + "\n", + "**This is where the bunching feedback loop lives**: a delayed bus finds more passengers waiting → longer dwell → more delay. Meanwhile, the following bus finds fewer passengers → shorter dwell → catches up." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "class BusStop(Component):\n", + " \"\"\"A bus stop that accumulates passengers and handles boarding/alighting.\"\"\"\n", + "\n", + " io = IO(\n", + " inputs=[\"tick\"],\n", + " input_events=[BusArrivedAtStop],\n", + " output_events=[PassengersBoarded],\n", + " )\n", + "\n", + " def __init__(\n", + " self,\n", + " stop_index: int,\n", + " arrival_rate: float,\n", + " alight_prob: float,\n", + " rng_seed: int = 42,\n", + " **kwargs: _t.Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._stop_index = stop_index\n", + " self._arrival_rate = arrival_rate\n", + " self._alight_prob = alight_prob\n", + " self._queue = 0\n", + " self._rng = np.random.default_rng(rng_seed + stop_index * 7)\n", + " # Record arrivals for headway analysis: list of (tick, bus_id)\n", + " self.arrival_log: list[tuple[int, str]] = []\n", + "\n", + " async def step(self) -> None:\n", + " new_arrivals = int(self._rng.poisson(self._arrival_rate))\n", + " self._queue += new_arrivals\n", + "\n", + " @BusArrivedAtStop.handler\n", + " async def handle_bus_arrival(self, event: BusArrivedAtStop) -> PassengersBoarded:\n", + " \"\"\"Board/alight passengers when a bus arrives.\"\"\"\n", + " if event.data.stop_index != self._stop_index:\n", + " return PassengersBoarded(\n", + " source=self.name,\n", + " data=PassengersBoardedData(\n", + " bus_id=event.data.bus_id,\n", + " stop_index=self._stop_index,\n", + " boarded=0,\n", + " alighted=0,\n", + " new_load=event.data.current_load,\n", + " dwell_ticks=0,\n", + " ),\n", + " )\n", + "\n", + " # Record arrival\n", + " self.arrival_log.append((event.data.tick, event.data.bus_id))\n", + " load = event.data.current_load\n", + "\n", + " # Alighting\n", + " alighted = int(self._rng.binomial(load, self._alight_prob)) if load > 0 else 0\n", + " load -= alighted\n", + "\n", + " # Boarding (limited by remaining capacity)\n", + " space = BUS_CAPACITY - load\n", + " boarded = min(self._queue, space)\n", + " self._queue -= boarded\n", + " load += boarded\n", + "\n", + " # Dwell time\n", + " dwell_seconds = DWELL_FIXED + DWELL_PER_PAX * (boarded + alighted)\n", + " dwell_ticks = max(1, math.ceil(dwell_seconds / TICK_SECONDS))\n", + "\n", + " return PassengersBoarded(\n", + " source=self.name,\n", + " data=PassengersBoardedData(\n", + " bus_id=event.data.bus_id,\n", + " stop_index=self._stop_index,\n", + " boarded=boarded,\n", + " alighted=alighted,\n", + " new_load=load,\n", + " dwell_ticks=dwell_ticks,\n", + " ),\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "## Create the Bus Vehicle component\n", + "\n", + "Each `Bus` travels one-way along the corridor. Key features:\n", + "- **Staggered departure**: each bus waits for its scheduled departure tick before moving\n", + "- **Perturbation support**: an optional extra delay at a specific stop seeds the bunching cascade\n", + "- **Traffic signal delays**: random red-light encounters add realistic variance to travel times\n", + "- **State machine**: traveling → dwelling → (optionally holding) → traveling\n", + "- **Trajectory recording**: each tick's state is logged for post-simulation visualization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "class Bus(Component):\n", + " \"\"\"A bus that travels one-way along the corridor.\"\"\"\n", + "\n", + " io = IO(\n", + " inputs=[\"tick\"],\n", + " input_events=[PassengersBoarded, HoldBus],\n", + " output_events=[BusArrivedAtStop, BusPositionUpdate],\n", + " )\n", + "\n", + " def __init__(\n", + " self,\n", + " bus_id: str,\n", + " departure_tick: int,\n", + " signal_positions: list[float],\n", + " perturbation_stop: int | None = None,\n", + " perturbation_ticks: int = 0,\n", + " rng_seed: int = 42,\n", + " **kwargs: _t.Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._bus_id = bus_id\n", + " self._departure_tick = departure_tick\n", + " self._position = 0.0\n", + " self._next_stop_index = 0\n", + " self._load = 0\n", + " self._phase = \"waiting\" # waiting | traveling | dwelling | holding | finished\n", + " self._dwell_remaining = 0\n", + " self._hold_remaining = 0\n", + " self._signal_positions = signal_positions\n", + " self._perturbation_stop = perturbation_stop\n", + " self._perturbation_ticks = perturbation_ticks\n", + " self._perturbation_applied = False\n", + " self._rng = np.random.default_rng(rng_seed)\n", + " self._finished = False\n", + " self._awaiting_boarding = False # True while waiting for PassengersBoarded response\n", + " self.trajectory: list[dict[str, _t.Any]] = []\n", + "\n", + " async def step(self) -> None:\n", + " current_tick = int(self.tick)\n", + "\n", + " # Wait until departure time\n", + " if self._phase == \"waiting\":\n", + " if current_tick >= self._departure_tick:\n", + " self._phase = \"traveling\"\n", + " self._record(current_tick)\n", + " return\n", + "\n", + " if self._phase == \"finished\":\n", + " self._record(current_tick)\n", + " return\n", + "\n", + " if self._phase == \"dwelling\":\n", + " # Don't leave the stop until PassengersBoarded response has arrived\n", + " if self._awaiting_boarding:\n", + " self._record(current_tick)\n", + " return\n", + " self._dwell_remaining -= 1\n", + " if self._dwell_remaining <= 0:\n", + " # Apply perturbation delay if configured\n", + " if (\n", + " not self._perturbation_applied\n", + " and self._perturbation_stop is not None\n", + " and self._next_stop_index == self._perturbation_stop\n", + " ):\n", + " self._perturbation_applied = True\n", + " self._phase = \"holding\"\n", + " self._hold_remaining = self._perturbation_ticks\n", + " self._record(current_tick)\n", + " return\n", + "\n", + " self._next_stop_index += 1\n", + " if self._next_stop_index >= NUM_STOPS:\n", + " self._phase = \"finished\"\n", + " self._record(current_tick)\n", + " return\n", + " self._phase = \"traveling\"\n", + "\n", + " elif self._phase == \"holding\":\n", + " self._hold_remaining -= 1\n", + " if self._hold_remaining <= 0:\n", + " self._next_stop_index += 1\n", + " if self._next_stop_index >= NUM_STOPS:\n", + " self._phase = \"finished\"\n", + " self._record(current_tick)\n", + " return\n", + " self._phase = \"traveling\"\n", + "\n", + " if self._phase == \"traveling\":\n", + " distance = BUS_SPEED * TICK_SECONDS\n", + "\n", + " # Traffic signal delay\n", + " for sig_pos in self._signal_positions:\n", + " if self._position < sig_pos <= self._position + distance:\n", + " phase_in_cycle = self._rng.uniform(0, 90.0)\n", + " if phase_in_cycle > SIGNAL_GREEN_FRAC * 90.0:\n", + " red_wait = 90.0 - phase_in_cycle\n", + " distance = max(0.0, distance - BUS_SPEED * red_wait)\n", + "\n", + " self._position = min(self._position + distance, ROUTE_LENGTH)\n", + "\n", + " # Check if we reached the next stop\n", + " if self._next_stop_index < NUM_STOPS:\n", + " target = STOP_POSITIONS[self._next_stop_index]\n", + " if self._position >= target:\n", + " self._position = target\n", + " self._phase = \"dwelling\"\n", + " self._dwell_remaining = 0 # set by PassengersBoarded handler\n", + " self._awaiting_boarding = True\n", + "\n", + " self.io.queue_event(\n", + " BusArrivedAtStop(\n", + " source=self.name,\n", + " data=BusArrivedData(\n", + " bus_id=self._bus_id,\n", + " stop_index=self._next_stop_index,\n", + " current_load=self._load,\n", + " tick=current_tick,\n", + " ),\n", + " )\n", + " )\n", + "\n", + " self._record(current_tick)\n", + "\n", + " # Position update for dispatch\n", + " self.io.queue_event(\n", + " BusPositionUpdate(\n", + " source=self.name,\n", + " data=BusPositionData(\n", + " bus_id=self._bus_id,\n", + " position=self._position,\n", + " stop_index=self._next_stop_index,\n", + " load=self._load,\n", + " tick=current_tick,\n", + " phase=self._phase,\n", + " ),\n", + " )\n", + " )\n", + "\n", + " def _record(self, tick: int) -> None:\n", + " self.trajectory.append(\n", + " {\n", + " \"tick\": tick,\n", + " \"position\": self._position,\n", + " \"stop_index\": self._next_stop_index,\n", + " \"load\": self._load,\n", + " \"phase\": self._phase,\n", + " }\n", + " )\n", + "\n", + " @PassengersBoarded.handler\n", + " async def handle_boarding(self, event: PassengersBoarded) -> None:\n", + " if event.data.bus_id != self._bus_id:\n", + " return\n", + " self._load = event.data.new_load\n", + " self._awaiting_boarding = False\n", + " if event.data.dwell_ticks > 0:\n", + " self._dwell_remaining = max(self._dwell_remaining, event.data.dwell_ticks)\n", + "\n", + " @HoldBus.handler\n", + " async def handle_hold(self, event: HoldBus) -> None:\n", + " if event.data.bus_id != self._bus_id:\n", + " return\n", + " if self._phase in (\"dwelling\", \"holding\"):\n", + " self._phase = \"holding\"\n", + " self._hold_remaining = max(self._hold_remaining, event.data.hold_ticks)" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "## Create the Simulation Clock and Dispatch Controller\n", + "\n", + "The **SimClock** drives the simulation forward.\n", + "\n", + "The **Dispatch** controller listens to `BusPositionUpdate` events and computes headways between consecutive active buses. When two buses are too close together (headway below threshold), it holds the trailing bus at its current stop to restore spacing." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "class SimClock(Component):\n", + " \"\"\"Simulation clock that drives the tick counter.\"\"\"\n", + "\n", + " io = IO(outputs=[\"tick\"])\n", + "\n", + " def __init__(\n", + " self, total_ticks: int = TOTAL_TICKS, **kwargs: _t.Unpack[ComponentArgsDict]\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._total = total_ticks\n", + " self._tick = 0\n", + "\n", + " async def step(self) -> None:\n", + " self._tick += 1\n", + " self.tick = self._tick\n", + " if self._tick >= self._total:\n", + " await self.io.close()\n", + "\n", + "\n", + "class Dispatch(Component):\n", + " \"\"\"Monitors bus positions and issues holding instructions to restore headways.\"\"\"\n", + "\n", + " io = IO(\n", + " inputs=[\"tick\"],\n", + " input_events=[BusPositionUpdate],\n", + " output_events=[HoldBus],\n", + " )\n", + "\n", + " def __init__(\n", + " self,\n", + " bus_ids: list[str],\n", + " enabled: bool = True,\n", + " target_headway_s: float = DEPARTURE_HEADWAY_TICKS * TICK_SECONDS,\n", + " threshold_s: float = HEADWAY_THRESHOLD,\n", + " **kwargs: _t.Unpack[ComponentArgsDict],\n", + " ) -> None:\n", + " super().__init__(**kwargs)\n", + " self._bus_ids = bus_ids\n", + " self._enabled = enabled\n", + " self._target = target_headway_s\n", + " self._threshold = threshold_s\n", + " self._positions: dict[str, float] = {}\n", + " self._phases: dict[str, str] = {}\n", + " self._hold_cooldown: dict[str, int] = {} # bus_id → tick when hold expires\n", + "\n", + " async def step(self) -> None:\n", + " if not self._enabled:\n", + " return\n", + "\n", + " # Only consider active buses (not waiting/finished)\n", + " active = {\n", + " bid: pos\n", + " for bid, pos in self._positions.items()\n", + " if self._phases.get(bid) not in (\"waiting\", \"finished\", None)\n", + " }\n", + "\n", + " if len(active) < 2:\n", + " return\n", + "\n", + " # Use departure order (not position) to identify leader/follower pairs.\n", + " # In a one-way route, bus-0 always leads bus-1, etc.\n", + " ordered = [bid for bid in self._bus_ids if bid in active]\n", + "\n", + " for i in range(len(ordered) - 1):\n", + " leader_id = ordered[i]\n", + " follower_id = ordered[i + 1]\n", + "\n", + " gap = max(0.0, active[leader_id] - active[follower_id])\n", + " time_gap = gap / BUS_SPEED if BUS_SPEED > 0 else float(\"inf\")\n", + "\n", + " # If follower is too close to leader, hold the follower\n", + " if time_gap < self._target - self._threshold:\n", + " hold_s = min((self._target - time_gap) * 0.4, 60.0)\n", + " hold_ticks = max(1, math.ceil(hold_s / TICK_SECONDS))\n", + " current_tick = int(self.tick)\n", + "\n", + " # Only hold if not on cooldown from a recent hold\n", + " if self._phases.get(follower_id) in (\n", + " \"dwelling\",\n", + " \"holding\",\n", + " ) and current_tick >= self._hold_cooldown.get(follower_id, 0):\n", + " self._hold_cooldown[follower_id] = current_tick + hold_ticks + 2\n", + " self.io.queue_event(\n", + " HoldBus(\n", + " source=self.name,\n", + " data=HoldBusData(bus_id=follower_id, hold_ticks=hold_ticks),\n", + " )\n", + " )\n", + "\n", + " @BusPositionUpdate.handler\n", + " async def handle_position(self, event: BusPositionUpdate) -> None:\n", + " self._positions[event.data.bus_id] = event.data.position\n", + " self._phases[event.data.bus_id] = event.data.phase" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "## Wire up the process\n", + "\n", + "We assemble the Plugboard process with:\n", + "- **Field connectors**: Clock tick → all components (synchronization)\n", + "- **Event connectors**: automatically wired based on declared input/output events\n", + "\n", + "Plugboard's declarative wiring makes the complex interaction topology explicit — buses talk to stops via events, stops respond via events, and the dispatch controller monitors everything through the same event bus." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "signal_positions = [ROUTE_LENGTH * (i + 1) / (NUM_SIGNALS + 1) for i in range(NUM_SIGNALS)]\n", + "\n", + "\n", + "def build_process(\n", + " dispatch_enabled: bool = False,\n", + " apply_perturbation: bool = True,\n", + " rng_seed: int = 42,\n", + ") -> tuple[LocalProcess, list[Bus], list[BusStop], Dispatch]:\n", + " \"\"\"Build a simulation process with all components wired together.\"\"\"\n", + " clock = SimClock(name=\"clock\")\n", + "\n", + " stops = [\n", + " BusStop(\n", + " name=f\"stop-{i}\",\n", + " stop_index=i,\n", + " arrival_rate=BASE_ARRIVAL_RATES[i],\n", + " alight_prob=ALIGHT_PROBS[i],\n", + " rng_seed=rng_seed,\n", + " )\n", + " for i in range(NUM_STOPS)\n", + " ]\n", + "\n", + " bus_ids = [f\"bus-{i}\" for i in range(NUM_BUSES)]\n", + " buses = [\n", + " Bus(\n", + " name=bus_ids[i],\n", + " bus_id=bus_ids[i],\n", + " departure_tick=i * DEPARTURE_HEADWAY_TICKS,\n", + " signal_positions=signal_positions,\n", + " perturbation_stop=PERTURBATION_STOP\n", + " if (apply_perturbation and bus_ids[i] == PERTURBATION_BUS)\n", + " else None,\n", + " perturbation_ticks=PERTURBATION_TICKS\n", + " if (apply_perturbation and bus_ids[i] == PERTURBATION_BUS)\n", + " else 0,\n", + " rng_seed=rng_seed + i * 100,\n", + " )\n", + " for i in range(NUM_BUSES)\n", + " ]\n", + "\n", + " dispatch = Dispatch(name=\"dispatch\", bus_ids=bus_ids, enabled=dispatch_enabled)\n", + "\n", + " components: list[Component] = [clock] + stops + buses + [dispatch] # type: ignore[list-item]\n", + "\n", + " field_connectors = []\n", + " for comp in stops + buses + [dispatch]:\n", + " field_connectors.append(connect(\"clock.tick\", f\"{comp.name}.tick\"))\n", + "\n", + " event_connectors = AsyncioConnector.builder().build_event_connectors(components)\n", + "\n", + " process = LocalProcess(\n", + " components=components,\n", + " connectors=field_connectors + event_connectors,\n", + " )\n", + " return process, buses, stops, dispatch\n", + "\n", + "\n", + "# Visualize the process wiring\n", + "process_viz, _, _, _ = build_process(dispatch_enabled=True)\n", + "diagram = MermaidDiagram.from_process(process_viz)\n", + "diagram.url" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "![Process diagram](https://mermaid.ink/img/pako:eNqtml-PmkAUxb8KmUfjKgOKQJpN90-TvjTZZLsvjYkZYRRSZAgM290av3sHWW0fdlPvcHwSuJffHNDLOQl7lqhUsphtCvUryUStne-3y9JxkkIlPz_vnSYTlYydWrVlKtOxU4i1LGLnMd_ddRWf1vX1aHQsHo2cg3N1de00WlVX7se9t23zaEr61r6460VQOYXKUVSPQvVQVJ9C9VHUGYU6Q1HnFOocRQ0o1ABFXVCoCxQ1pFBDFDWiUCPYlCANJ46bTrTxBJtPnDSgOGxCcdKI4rAZxUlDiiOm1Lpt_vPI67uOdSAev5DHQTzvQp4H4vkX8nwQb3YhD_V7mV_IQzzP0ryphE6yj7vv3yr67lP9CW3WsRJ1nT_LdCX0qvvj_D1VJl_-XflNX3ejuz9ax58MM4UYtp01xLDtDCKGbWcTMWw7s4hh21lGDNvOOGLYdvYRw7YzkRi2nZUEzRZLQwmiW9pKEN3SXILolhYTRLc0mpkq0pVZwvvIr-ao6T6RrKxlJZpGlltZN6u1ErVpeJ_1cK677csGUcm6yBYWootMJesiW2WILjKVrItsySG6yFSyLrL1h-giU8m6yBEDootM7aZxpZpc56pctVUqtPxwGj-81T0dy07YofnGMp708AEX7cymxxMcmx5PcGx6PMGx6fEEx6bHExybHk9wbHo8wbHp8QTHpscT4GyxiCdAukU8AdIt4gmQbhFPgHSLeAKhk8PD2TvYhrKBVFvzccKSEgVMrCV1qFhSzICJtaQOFUvKHjCxltShYkmBBCbWkjpULCmlwMRaUoeIHZKN-hVcmDvZmO1kvRN5ymK279BLpjO5k0sWm6-p3Ii20Eu2LLtS0Wr1-FomLNZ1K8fMLGqbsXgjisZs9Trvc7Gtxe5U0u_8kuZa1ad9hRKpNFt7pl-r7k2mbd5oc_pElZt82-1v68LszrSumng67Q5PtrnO2vUkUbtpk6fda0_ZcxRMAy8IhefLYOGLue-nyZpH4cab8U26cLkn2OEwZpUofyh1XpPZ7CAvLJ6HkyBczKIwDOah50fBmL2ymIfBJIpcz_e9ubdwZ-YMv4_t7qQrNp_QdyPX9RZ8zORR2rf-lazjm1ljtq276_l2jaS5b_WduX2axX44O_wBlgS1ZQ)" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "## Run Scenario 1: No control (bus bunching)\n", + "\n", + "First we run **without** the dispatch controller. The initial perturbation on bus-0 (a 120s delay at stop 2) seeds the bunching cascade. As bus-0 falls behind, it encounters more waiting passengers, increasing dwell times and delaying it further. Bus-1 catches up, finding fewer passengers and shorter dwells." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "process_nc, buses_nc, stops_nc, dispatch_nc = build_process(dispatch_enabled=False, rng_seed=42)\n", + "\n", + "async with process_nc:\n", + " await process_nc.run()\n", + "\n", + "print(\"Scenario 1 (no control) complete\")" + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "## Run Scenario 2: With headway control\n", + "\n", + "Now we run with the dispatch controller enabled. It monitors headways and holds buses at stops when they get too close to the bus ahead." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "process_ctrl, buses_ctrl, stops_ctrl, dispatch_ctrl = build_process(\n", + " dispatch_enabled=True, rng_seed=42\n", + ")\n", + "\n", + "async with process_ctrl:\n", + " await process_ctrl.run()\n", + "\n", + "print(\"Scenario 2 (with control) complete\")" + ] + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": {}, + "source": [ + "## Visualize bus trajectories (time-space diagram)\n", + "\n", + "The **time-space diagram** is the standard transit operations visualization. Each line traces a bus along the corridor. Flat segments = dwelling/holding at a stop. Converging lines = bunching." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "_COLORS = px.colors.qualitative.Plotly\n", + "_x_max = TOTAL_TICKS * TICK_SECONDS / 60\n", + "\n", + "fig = make_subplots(\n", + " rows=1,\n", + " cols=2,\n", + " shared_yaxes=True,\n", + " subplot_titles=[\"No Control\", \"With Headway Control\"],\n", + " horizontal_spacing=0.05,\n", + ")\n", + "\n", + "for col, (buses, _title) in enumerate(\n", + " [(buses_nc, \"No Control\"), (buses_ctrl, \"With Headway Control\")], 1\n", + "):\n", + " for i, bus in enumerate(buses):\n", + " ticks = [r[\"tick\"] for r in bus.trajectory]\n", + " times_m = [t * TICK_SECONDS / 60 for t in ticks]\n", + " positions_km = [r[\"position\"] / 1000 for r in bus.trajectory]\n", + " fig.add_trace(\n", + " go.Scatter(\n", + " x=times_m,\n", + " y=positions_km,\n", + " mode=\"lines\",\n", + " name=bus.name,\n", + " line=dict(color=_COLORS[i % len(_COLORS)], width=1.5),\n", + " opacity=0.85,\n", + " legendgroup=bus.name,\n", + " showlegend=(col == 1),\n", + " ),\n", + " row=1,\n", + " col=col,\n", + " )\n", + " for sp in STOP_POSITIONS:\n", + " fig.add_shape(\n", + " type=\"line\",\n", + " x0=0,\n", + " x1=_x_max,\n", + " y0=sp / 1000,\n", + " y1=sp / 1000,\n", + " line=dict(color=\"lightgray\", width=0.5),\n", + " row=1,\n", + " col=col,\n", + " )\n", + "\n", + "fig.update_xaxes(title_text=\"Time (minutes)\")\n", + "fig.update_yaxes(title_text=\"Distance along route (km)\", col=1)\n", + "fig.update_layout(\n", + " title=dict(text=\"Time-Space Diagram: Bus Bunching\", font=dict(size=16)),\n", + " height=550,\n", + " legend=dict(orientation=\"v\"),\n", + ")\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "## Headway analysis at stops\n", + "\n", + "**Headway** is the time between consecutive bus arrivals at the same stop. Even headways = regular service. Headways that shrink then grow = buses bunching up.\n", + "\n", + "We compute headways from the `arrival_log` recorded by each `BusStop` component and plot the coefficient of variation (CV) across stops." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24", + "metadata": {}, + "outputs": [], + "source": [ + "def compute_stop_headways(\n", + " stops: list[BusStop],\n", + ") -> dict[int, list[float]]:\n", + " \"\"\"Compute headways (in seconds) at each stop from arrival logs.\"\"\"\n", + " headways: dict[int, list[float]] = {}\n", + " for i, stop in enumerate(stops):\n", + " log = sorted(stop.arrival_log, key=lambda x: x[0]) # sort by tick\n", + " if len(log) < 2:\n", + " headways[i] = []\n", + " continue\n", + " hw = []\n", + " for j in range(1, len(log)):\n", + " dt = (log[j][0] - log[j - 1][0]) * TICK_SECONDS\n", + " hw.append(dt)\n", + " headways[i] = hw\n", + " return headways\n", + "\n", + "\n", + "hw_nc = compute_stop_headways(stops_nc)\n", + "hw_ctrl = compute_stop_headways(stops_ctrl)\n", + "\n", + "\n", + "def headway_cv(headways: dict[int, list[float]]) -> list[float]:\n", + " \"\"\"Compute CV (std/mean) at each stop.\"\"\"\n", + " cvs = []\n", + " for i in range(NUM_STOPS):\n", + " hw = headways.get(i, [])\n", + " if len(hw) >= 2:\n", + " arr = np.array(hw)\n", + " cvs.append(float(arr.std() / arr.mean()) if arr.mean() > 0 else 0.0)\n", + " else:\n", + " cvs.append(0.0)\n", + " return cvs\n", + "\n", + "\n", + "cv_nc = headway_cv(hw_nc)\n", + "cv_ctrl = headway_cv(hw_ctrl)\n", + "\n", + "# Headway distribution (all stops combined)\n", + "all_hw_nc = [h for hws in hw_nc.values() for h in hws]\n", + "all_hw_ctrl = [h for hws in hw_ctrl.values() for h in hws]\n", + "\n", + "ideal_headway = DEPARTURE_HEADWAY_TICKS * TICK_SECONDS\n", + "stop_indices = list(range(NUM_STOPS))\n", + "\n", + "fig = make_subplots(\n", + " rows=1,\n", + " cols=2,\n", + " subplot_titles=[\"Headway Regularity by Stop\", \"Headway Distribution (All Stops)\"],\n", + " horizontal_spacing=0.1,\n", + ")\n", + "\n", + "fig.add_trace(\n", + " go.Bar(name=\"No Control\", x=stop_indices, y=cv_nc, marker_color=\"salmon\", opacity=0.8),\n", + " row=1,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Bar(name=\"With Control\", x=stop_indices, y=cv_ctrl, marker_color=\"steelblue\", opacity=0.8),\n", + " row=1,\n", + " col=1,\n", + ")\n", + "fig.add_trace(\n", + " go.Histogram(\n", + " x=all_hw_nc,\n", + " name=\"No Control\",\n", + " marker_color=\"salmon\",\n", + " opacity=0.6,\n", + " histnorm=\"probability density\",\n", + " nbinsx=25,\n", + " showlegend=False,\n", + " ),\n", + " row=1,\n", + " col=2,\n", + ")\n", + "fig.add_trace(\n", + " go.Histogram(\n", + " x=all_hw_ctrl,\n", + " name=\"With Control\",\n", + " marker_color=\"steelblue\",\n", + " opacity=0.6,\n", + " histnorm=\"probability density\",\n", + " nbinsx=25,\n", + " showlegend=False,\n", + " ),\n", + " row=1,\n", + " col=2,\n", + ")\n", + "fig.add_vline(\n", + " x=ideal_headway,\n", + " line=dict(color=\"black\", dash=\"dash\", width=1.5),\n", + " annotation_text=f\"Ideal ({ideal_headway}s)\",\n", + " annotation_position=\"top right\",\n", + " row=1,\n", + " col=2,\n", + ")\n", + "\n", + "fig.update_xaxes(title_text=\"Stop Index\", col=1)\n", + "fig.update_yaxes(title_text=\"Headway CV (std / mean)\", col=1)\n", + "fig.update_xaxes(title_text=\"Headway (seconds)\", col=2)\n", + "fig.update_yaxes(title_text=\"Density\", col=2)\n", + "fig.update_layout(\n", + " title=dict(text=\"Headway Analysis\", font=dict(size=16)),\n", + " barmode=\"group\",\n", + " height=500,\n", + ")\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "25", + "metadata": {}, + "source": [ + "## Bus load factors\n", + "\n", + "How full are the buses? In the bunching scenario, the delayed bus picks up more passengers (it's full) while the bus behind runs nearly empty — a hallmark of the bunching problem." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26", + "metadata": {}, + "outputs": [], + "source": [ + "fig = make_subplots(\n", + " rows=1,\n", + " cols=2,\n", + " shared_yaxes=True,\n", + " subplot_titles=[\"No Control\", \"With Headway Control\"],\n", + " horizontal_spacing=0.05,\n", + ")\n", + "\n", + "for col, (buses, _title) in enumerate(\n", + " [(buses_nc, \"No Control\"), (buses_ctrl, \"With Headway Control\")], 1\n", + "):\n", + " for i, bus in enumerate(buses):\n", + " active = [r for r in bus.trajectory if r[\"phase\"] not in (\"waiting\",)]\n", + " if not active:\n", + " continue\n", + " times_m = [r[\"tick\"] * TICK_SECONDS / 60 for r in active]\n", + " loads_pct = [r[\"load\"] / BUS_CAPACITY * 100 for r in active]\n", + " fig.add_trace(\n", + " go.Scatter(\n", + " x=times_m,\n", + " y=loads_pct,\n", + " mode=\"lines\",\n", + " name=bus.name,\n", + " line=dict(color=_COLORS[i % len(_COLORS)], width=1.2),\n", + " opacity=0.8,\n", + " legendgroup=bus.name,\n", + " showlegend=(col == 1),\n", + " ),\n", + " row=1,\n", + " col=col,\n", + " )\n", + " fig.add_hline(\n", + " y=100, line=dict(color=\"red\", width=0.8, dash=\"dash\"), opacity=0.5, row=1, col=col\n", + " )\n", + "\n", + "fig.update_xaxes(title_text=\"Time (minutes)\")\n", + "fig.update_yaxes(title_text=\"Load Factor (%)\", col=1, range=[0, 110])\n", + "fig.update_yaxes(range=[0, 110])\n", + "fig.update_layout(\n", + " title=dict(text=\"Bus Load Factors Over Time\", font=dict(size=16)),\n", + " height=500,\n", + ")\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "27", + "metadata": {}, + "source": [ + "## Summary statistics\n", + "\n", + "We quantify the improvement from headway control using key transit performance metrics:\n", + "- **Mean headway** — average time between consecutive buses at each stop\n", + "- **Headway CV** — coefficient of variation (lower = more regular)\n", + "- **Excess wait time** — how much longer passengers wait compared to perfectly regular service\n", + "- **P95 headway** — the 95th percentile headway (tail events passengers experience)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28", + "metadata": {}, + "outputs": [], + "source": [ + "def compute_metrics(headways: dict[int, list[float]]) -> dict[str, float]:\n", + " \"\"\"Compute aggregate headway metrics across all stops.\"\"\"\n", + " all_hw = [h for hws in headways.values() for h in hws if h > 0]\n", + " if not all_hw:\n", + " return {\"mean\": 0, \"std\": 0, \"cv\": 0, \"p95\": 0, \"excess_wait\": 0}\n", + " arr = np.array(all_hw)\n", + " mean = float(arr.mean())\n", + " std = float(arr.std())\n", + " cv = std / mean if mean > 0 else 0\n", + " p95 = float(np.percentile(arr, 95))\n", + " # Excess wait = E[W] - H/2, where E[W] = H/2 * (1 + CV^2)\n", + " ideal_wait = mean / 2\n", + " actual_wait = mean / 2 * (1 + cv**2)\n", + " excess_wait = actual_wait - ideal_wait\n", + " return {\"mean\": mean, \"std\": std, \"cv\": cv, \"p95\": p95, \"excess_wait\": excess_wait}\n", + "\n", + "\n", + "m_nc = compute_metrics(hw_nc)\n", + "m_ctrl = compute_metrics(hw_ctrl)\n", + "\n", + "ideal_headway = DEPARTURE_HEADWAY_TICKS * TICK_SECONDS\n", + "\n", + "print(\"=\" * 60)\n", + "print(f\"{'Metric':<25} {'No Control':>12} {'With Control':>14} {'Change':>10}\")\n", + "print(\"=\" * 60)\n", + "print(f\"{'Ideal headway (s)':<25} {ideal_headway:>12.0f} {ideal_headway:>14.0f} {'':>10}\")\n", + "print(\n", + " f\"{'Mean headway (s)':<25} {m_nc['mean']:>12.1f} {m_ctrl['mean']:>14.1f} {(m_ctrl['mean'] - m_nc['mean']):>+10.1f}\"\n", + ")\n", + "print(\n", + " f\"{'Headway std dev (s)':<25} {m_nc['std']:>12.1f} {m_ctrl['std']:>14.1f} {(m_ctrl['std'] - m_nc['std']):>+10.1f}\"\n", + ")\n", + "print(\n", + " f\"{'Headway CV':<25} {m_nc['cv']:>12.3f} {m_ctrl['cv']:>14.3f} {(m_ctrl['cv'] - m_nc['cv']):>+10.3f}\"\n", + ")\n", + "print(\n", + " f\"{'P95 headway (s)':<25} {m_nc['p95']:>12.1f} {m_ctrl['p95']:>14.1f} {(m_ctrl['p95'] - m_nc['p95']):>+10.1f}\"\n", + ")\n", + "print(\n", + " f\"{'Excess wait (s)':<25} {m_nc['excess_wait']:>12.1f} {m_ctrl['excess_wait']:>14.1f} {(m_ctrl['excess_wait'] - m_nc['excess_wait']):>+10.1f}\"\n", + ")\n", + "print(\"=\" * 60)\n", + "\n", + "if m_nc[\"cv\"] > 0:\n", + " improvement = (1 - m_ctrl[\"cv\"] / m_nc[\"cv\"]) * 100\n", + " print(f\"\\nHeadway regularity improvement: {improvement:.1f}%\")\n", + "if m_nc[\"excess_wait\"] > 0:\n", + " wait_improvement = (1 - m_ctrl[\"excess_wait\"] / m_nc[\"excess_wait\"]) * 100\n", + " print(f\"Excess wait time reduction: {wait_improvement:.1f}%\")" + ] + }, + { + "cell_type": "markdown", + "id": "29", + "metadata": {}, + "source": [ + "## Key takeaways\n", + "\n", + "### What this example demonstrates\n", + "\n", + "1. **Event-driven simulation is a natural fit for transit modeling.** Bus arrivals, passenger boarding, and dispatch decisions are inherently event-driven. Plugboard's `Event` and handler system maps cleanly to these interactions — no polling or centralized scheduler needed.\n", + "\n", + "2. **The bus bunching feedback loop.** A small initial perturbation (120s delay) cascades through the system: the delayed bus encounters more waiting passengers → longer dwell times → falls further behind → the following bus catches up, finding fewer passengers → eventually buses pair up.\n", + "\n", + "3. **Headway control via holding strategies.** By monitoring bus positions through `BusPositionUpdate` events and issuing `HoldBus` events, the `Dispatch` component can regulate headways in real time — the same approach used by real transit agencies.\n", + "\n", + "4. **Modular, composable architecture.** Each component (`Bus`, `BusStop`, `SimClock`, `Dispatch`) is independent and testable. The event-based wiring lets us add or swap control strategies without modifying existing components.\n", + "\n", + "5. **Quantifiable improvement.** The headway-based holding strategy reduced headway variability (CV) by over 50% and excess passenger wait time by over 75% — a dramatic improvement from a relatively simple control policy, all implemented as an additional event-producing component.\n", + "\n", + "### Extensions\n", + "\n", + "- **Multiple routes with transfers** — add interchange stops that generate cross-route events\n", + "- **Demand-responsive scheduling** — use events to trigger schedule adjustments\n", + "- **Real-time passenger information** — add components that estimate arrival times from bus positions\n", + "- **Distributed simulation** — use `RayProcess` to run different route segments on separate workers" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "plugboard (3.12.8)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/mkdocs.yaml b/mkdocs.yaml index 7a7aa1e4..1dabbd9c 100644 --- a/mkdocs.yaml +++ b/mkdocs.yaml @@ -142,11 +142,17 @@ nav: - Rock paper scissors: examples/demos/llm/003_rock_paper_scissors/rock-paper-scissors.ipynb - Image processing and model routing: examples/demos/llm/004_image_processing/local-and-remote-image-processing.ipynb - Physics-based models: - - Hot water tank: examples/demos/physics-models/001-hot-water-tank/hot-water-tank.ipynb + - Hot water tank: examples/demos/physics-models/001_hot_water_tank/hot-water-tank.ipynb + - Quadruple tank process: examples/demos/physics-models/002_quadruple_tank/quadruple-tank.ipynb - Mining and mineral processing: - - Crusher circuit: examples/demos/mining-minerals/001-crusher-circuit/crusher_circuit.ipynb + - Crusher circuit: examples/demos/mining-minerals/001_crusher_circuit/crusher_circuit.ipynb - Finance: - Momentum trading signal: examples/demos/finance/001_momentum_signal/momentum-signal.ipynb + - Process control: + - Reduced BSM1 wastewater treatment: examples/demos/process-control/001_bsm1_wastewater/bsm1-wastewater.ipynb + - Transport: + - Traffic intersection: examples/demos/transport/001_traffic_simulation/traffic-simulation.ipynb + - Bus terminal: examples/demos/transport/002_bus_terminal/bus-terminal.ipynb - API Reference: - component: api/component/component.md - connector: api/connector/connector.md