Quadruple Tank Process¶

This notebook implements the Quadruple Tank Process — a classic MIMO (multi-input multi-output) benchmark problem from control theory, originally described by Johansson (2000).

The system consists of four interconnected water tanks arranged in two columns:

Lower tanks (Tank 1 & Tank 2): the controlled outputs
Upper tanks (Tank 3 & Tank 4): fed by pumps, drain into the lower tanks in a cross-coupled fashion
Two pumps with three-way valves split flow between upper and lower tanks

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.

We'll build the nonlinear process model in Plugboard, add decentralized PI controllers, run closed-loop simulations, and visualize the results.

Install and import dependencies¶

In [ ]:

Copied!

# Install plugboard and dependencies for Google Colab
!pip install -q plugboard plotly
# Install plugboard and dependencies for Google Colab
!pip install -q plugboard plotly

In [ ]:

Copied!





import math
import typing as _t

import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots

from plugboard.component import Component
from plugboard.component import IOController as IO
from plugboard.connector import AsyncioConnector
from plugboard.library import FileWriter
from plugboard.process import LocalProcess
from plugboard.schemas import ComponentArgsDict, ConnectorSpec
import math
import typing as _t

import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots

from plugboard.component import Component
from plugboard.component import IOController as IO
from plugboard.connector import AsyncioConnector
from plugboard.library import FileWriter
from plugboard.process import LocalProcess
from plugboard.schemas import ComponentArgsDict, ConnectorSpec

Define quadruple tank process parameters¶

We use the nominal parameter values from Johansson (2000). The key parameters are:

Parameter	Description	Value	Unit
$A_1, A_3$	Cross-section of Tanks 1, 3	28	cm²
$A_2, A_4$	Cross-section of Tanks 2, 4	32	cm²
$a_1, a_3$	Outlet hole area, Tanks 1, 3	0.071	cm²
$a_2, a_4$	Outlet hole area, Tanks 2, 4	0.057	cm²
$k_1$	Pump 1 flow constant	3.33	cm³/Vs
$k_2$	Pump 2 flow constant	3.35	cm³/Vs
$g$	Gravity	981	cm/s²
$\gamma_1, \gamma_2$	Valve split ratios	varies	—

When $\gamma_1 + \gamma_2 > 1$, the system is minimum phase (easier to control). When $\gamma_1 + \gamma_2 < 1$, the system is non-minimum phase (harder to control, exhibits inverse response).

In [ ]:

Copied!





# Physical parameters (Johansson 2000)
A1 = 28.0  # cm² - cross-section of Tank 1
A2 = 32.0  # cm² - cross-section of Tank 2
A3 = 28.0  # cm² - cross-section of Tank 3
A4 = 32.0  # cm² - cross-section of Tank 4

a1 = 0.071  # cm² - outlet hole area, Tank 1
a2 = 0.057  # cm² - outlet hole area, Tank 2
a3 = 0.071  # cm² - outlet hole area, Tank 3
a4 = 0.057  # cm² - outlet hole area, Tank 4

g = 981.0  # cm/s² - gravitational acceleration

k1 = 3.33  # cm³/Vs - pump 1 flow constant
k2 = 3.35  # cm³/Vs - pump 2 flow constant

# Minimum-phase valve settings (gamma1 + gamma2 > 1)
GAMMA1_MP = 0.70
GAMMA2_MP = 0.60

# Non-minimum-phase valve settings (gamma1 + gamma2 < 1)
GAMMA1_NMP = 0.43
GAMMA2_NMP = 0.34

# Nominal operating point pump voltages
V1_0 = 3.0  # V
V2_0 = 3.0  # V

# Steady-state levels for minimum-phase case (approximate, from Johansson 2000 Table I)
H1_0 = 12.4  # cm
H2_0 = 12.7  # cm
H3_0 = 1.8  # cm
H4_0 = 1.4  # cm

# Steady-state levels for non-minimum-phase case.
# Upper tank steady state from Torricelli balance: h_i = ((1-γ)*k*V / (a*√(2g)))²
# Lower tank steady state includes drainage from the upper tank.
H3_NMP = ((1.0 - GAMMA2_NMP) * k2 * V2_0 / (a3 * math.sqrt(2.0 * g))) ** 2  # ≈ 4.45 cm
H4_NMP = ((1.0 - GAMMA1_NMP) * k1 * V1_0 / (a4 * math.sqrt(2.0 * g))) ** 2  # ≈ 5.09 cm
H1_NMP = (
    (GAMMA1_NMP * k1 * V1_0 + a3 * math.sqrt(2.0 * g * H3_NMP)) / (a1 * math.sqrt(2.0 * g))
) ** 2
H2_NMP = (
    (GAMMA2_NMP * k2 * V2_0 + a4 * math.sqrt(2.0 * g * H4_NMP)) / (a2 * math.sqrt(2.0 * g))
) ** 2

# Simulation parameters
DT = 1.0  # s - time step
T_TOTAL = 600  # s - total simulation time
N_STEPS = int(T_TOTAL / DT)
# Physical parameters (Johansson 2000)
A1 = 28.0  # cm² - cross-section of Tank 1
A2 = 32.0  # cm² - cross-section of Tank 2
A3 = 28.0  # cm² - cross-section of Tank 3
A4 = 32.0  # cm² - cross-section of Tank 4

a1 = 0.071  # cm² - outlet hole area, Tank 1
a2 = 0.057  # cm² - outlet hole area, Tank 2
a3 = 0.071  # cm² - outlet hole area, Tank 3
a4 = 0.057  # cm² - outlet hole area, Tank 4

g = 981.0  # cm/s² - gravitational acceleration

k1 = 3.33  # cm³/Vs - pump 1 flow constant
k2 = 3.35  # cm³/Vs - pump 2 flow constant

# Minimum-phase valve settings (gamma1 + gamma2 > 1)
GAMMA1_MP = 0.70
GAMMA2_MP = 0.60

# Non-minimum-phase valve settings (gamma1 + gamma2 < 1)
GAMMA1_NMP = 0.43
GAMMA2_NMP = 0.34

# Nominal operating point pump voltages
V1_0 = 3.0  # V
V2_0 = 3.0  # V

# Steady-state levels for minimum-phase case (approximate, from Johansson 2000 Table I)
H1_0 = 12.4  # cm
H2_0 = 12.7  # cm
H3_0 = 1.8  # cm
H4_0 = 1.4  # cm

# Steady-state levels for non-minimum-phase case.
# Upper tank steady state from Torricelli balance: h_i = ((1-γ)*k*V / (a*√(2g)))²
# Lower tank steady state includes drainage from the upper tank.
H3_NMP = ((1.0 - GAMMA2_NMP) * k2 * V2_0 / (a3 * math.sqrt(2.0 * g))) ** 2  # ≈ 4.45 cm
H4_NMP = ((1.0 - GAMMA1_NMP) * k1 * V1_0 / (a4 * math.sqrt(2.0 * g))) ** 2  # ≈ 5.09 cm
H1_NMP = (
    (GAMMA1_NMP * k1 * V1_0 + a3 * math.sqrt(2.0 * g * H3_NMP)) / (a1 * math.sqrt(2.0 * g))
) ** 2
H2_NMP = (
    (GAMMA2_NMP * k2 * V2_0 + a4 * math.sqrt(2.0 * g * H4_NMP)) / (a2 * math.sqrt(2.0 * g))
) ** 2

# Simulation parameters
DT = 1.0  # s - time step
T_TOTAL = 600  # s - total simulation time
N_STEPS = int(T_TOTAL / DT)

Build the tank components¶

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.

Lower tanks (Tank 1, Tank 2) receive direct pump flow and drainage from the upper tanks:

$$\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$$

Upper tanks (Tank 3, Tank 4) receive only the cross-coupled pump flow:

$$\frac{dh_3}{dt} = -\frac{a_3}{A_3}\sqrt{2g\,h_3} + \frac{(1-\gamma_2) k_2}{A_3}\,v_2$$

We use a simple Euler integration step with time step dt.

In [ ]:

Copied!





class LowerTank(Component):
    """Lower tank: receives direct pump flow and drainage from an upper tank."""

    io = IO(inputs=["prev_h", "upper_h", "pump_voltage"], outputs=["h"])

    def __init__(
        self,
        tank_area: float = A1,
        outlet_area: float = a1,
        upper_outlet_area: float = a3,
        pump_gain: float = k1,
        gamma: float = GAMMA1_MP,
        gravity: float = g,
        dt: float = DT,
        **kwargs: _t.Unpack[ComponentArgsDict],
    ) -> None:
        super().__init__(**kwargs)
        self._A = tank_area
        self._a = outlet_area
        self._a_upper = upper_outlet_area
        self._k = pump_gain
        self._gamma = gamma
        self._dt = dt
        self._sqrt_2g = math.sqrt(2.0 * gravity)

    async def step(self) -> None:
        outflow = self._a * self._sqrt_2g * math.sqrt(max(self.prev_h, 0.0))
        inflow_upper = self._a_upper * self._sqrt_2g * math.sqrt(max(self.upper_h, 0.0))
        inflow_pump = self._gamma * self._k * self.pump_voltage

        dh = (-outflow + inflow_upper + inflow_pump) / self._A
        self.h = max(self.prev_h + dh * self._dt, 0.0)


class UpperTank(Component):
    """Upper tank: receives cross-coupled pump flow only."""

    io = IO(inputs=["prev_h", "pump_voltage"], outputs=["h"])

    def __init__(
        self,
        tank_area: float = A3,
        outlet_area: float = a3,
        pump_gain: float = k2,
        gamma_complement: float = 1.0 - GAMMA2_MP,
        gravity: float = g,
        dt: float = DT,
        **kwargs: _t.Unpack[ComponentArgsDict],
    ) -> None:
        super().__init__(**kwargs)
        self._A = tank_area
        self._a = outlet_area
        self._k = pump_gain
        self._gamma_c = gamma_complement
        self._dt = dt
        self._sqrt_2g = math.sqrt(2.0 * gravity)

    async def step(self) -> None:
        outflow = self._a * self._sqrt_2g * math.sqrt(max(self.prev_h, 0.0))
        inflow_pump = self._gamma_c * self._k * self.pump_voltage

        dh = (-outflow + inflow_pump) / self._A
        self.h = max(self.prev_h + dh * self._dt, 0.0)
class LowerTank(Component):
    """Lower tank: receives direct pump flow and drainage from an upper tank."""

    io = IO(inputs=["prev_h", "upper_h", "pump_voltage"], outputs=["h"])

    def __init__(
        self,
        tank_area: float = A1,
        outlet_area: float = a1,
        upper_outlet_area: float = a3,
        pump_gain: float = k1,
        gamma: float = GAMMA1_MP,
        gravity: float = g,
        dt: float = DT,
        **kwargs: _t.Unpack[ComponentArgsDict],
    ) -> None:
        super().__init__(**kwargs)
        self._A = tank_area
        self._a = outlet_area
        self._a_upper = upper_outlet_area
        self._k = pump_gain
        self._gamma = gamma
        self._dt = dt
        self._sqrt_2g = math.sqrt(2.0 * gravity)

    async def step(self) -> None:
        outflow = self._a * self._sqrt_2g * math.sqrt(max(self.prev_h, 0.0))
        inflow_upper = self._a_upper * self._sqrt_2g * math.sqrt(max(self.upper_h, 0.0))
        inflow_pump = self._gamma * self._k * self.pump_voltage

        dh = (-outflow + inflow_upper + inflow_pump) / self._A
        self.h = max(self.prev_h + dh * self._dt, 0.0)


class UpperTank(Component):
    """Upper tank: receives cross-coupled pump flow only."""

    io = IO(inputs=["prev_h", "pump_voltage"], outputs=["h"])

    def __init__(
        self,
        tank_area: float = A3,
        outlet_area: float = a3,
        pump_gain: float = k2,
        gamma_complement: float = 1.0 - GAMMA2_MP,
        gravity: float = g,
        dt: float = DT,
        **kwargs: _t.Unpack[ComponentArgsDict],
    ) -> None:
        super().__init__(**kwargs)
        self._A = tank_area
        self._a = outlet_area
        self._k = pump_gain
        self._gamma_c = gamma_complement
        self._dt = dt
        self._sqrt_2g = math.sqrt(2.0 * gravity)

    async def step(self) -> None:
        outflow = self._a * self._sqrt_2g * math.sqrt(max(self.prev_h, 0.0))
        inflow_pump = self._gamma_c * self._k * self.pump_voltage

        dh = (-outflow + inflow_pump) / self._A
        self.h = max(self.prev_h + dh * self._dt, 0.0)

Build the PI Controller component¶

We use decentralized PI control — two independent loops:

Loop 1: $v_1$ controls $h_1$
Loop 2: $v_2$ controls $h_2$

Each controller computes: $$u(t) = v_0 + K_p \cdot e(t) + K_i \int_0^t e(\tau)\,d\tau$$

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.

In [ ]:

Copied!





class PIController(Component):
    """Decentralized PI controller with anti-windup for two tank levels."""

    io = IO(
        inputs=["h1", "h2", "sp_h1", "sp_h2"],
        outputs=["v1", "v2"],
    )

    def __init__(
        self,
        kp: tuple[float, float] = (3.0, 2.7),
        ki: tuple[float, float] = (0.1, 0.068),
        v_min: float = 0.0,
        v_max: float = 10.0,
        v0: tuple[float, float] = (V1_0, V2_0),
        dt: float = DT,
        **kwargs: _t.Unpack[ComponentArgsDict],
    ) -> None:
        super().__init__(**kwargs)
        self._kp = kp
        self._ki = ki
        self._v_min = v_min
        self._v_max = v_max
        self._v0 = v0
        self._dt = dt
        self._integral = [0.0, 0.0]

    async def step(self) -> None:
        errors = [self.sp_h1 - self.h1, self.sp_h2 - self.h2]

        outputs = []
        for i, (kp, ki, v0, e) in enumerate(zip(self._kp, self._ki, self._v0, errors)):
            u_raw = v0 + kp * e + ki * self._integral[i]
            u_clamped = max(self._v_min, min(self._v_max, u_raw))
            # Anti-windup: only integrate when the output is not saturated
            if u_raw == u_clamped:
                self._integral[i] += e * self._dt
            outputs.append(u_clamped)

        # PI output = feedforward (steady-state voltage) + feedback
        self.v1, self.v2 = outputs[0], outputs[1]
class PIController(Component):
    """Decentralized PI controller with anti-windup for two tank levels."""

    io = IO(
        inputs=["h1", "h2", "sp_h1", "sp_h2"],
        outputs=["v1", "v2"],
    )

    def __init__(
        self,
        kp: tuple[float, float] = (3.0, 2.7),
        ki: tuple[float, float] = (0.1, 0.068),
        v_min: float = 0.0,
        v_max: float = 10.0,
        v0: tuple[float, float] = (V1_0, V2_0),
        dt: float = DT,
        **kwargs: _t.Unpack[ComponentArgsDict],
    ) -> None:
        super().__init__(**kwargs)
        self._kp = kp
        self._ki = ki
        self._v_min = v_min
        self._v_max = v_max
        self._v0 = v0
        self._dt = dt
        self._integral = [0.0, 0.0]

    async def step(self) -> None:
        errors = [self.sp_h1 - self.h1, self.sp_h2 - self.h2]

        outputs = []
        for i, (kp, ki, v0, e) in enumerate(zip(self._kp, self._ki, self._v0, errors)):
            u_raw = v0 + kp * e + ki * self._integral[i]
            u_clamped = max(self._v_min, min(self._v_max, u_raw))
            # Anti-windup: only integrate when the output is not saturated
            if u_raw == u_clamped:
                self._integral[i] += e * self._dt
            outputs.append(u_clamped)

        # PI output = feedforward (steady-state voltage) + feedback
        self.v1, self.v2 = outputs[0], outputs[1]

Build the Setpoint Generator component¶

The setpoint generator produces step changes:

$t=0$: both setpoints at steady-state levels ($h_1^0 = 12.4$ cm, $h_2^0 = 12.7$ cm)
$t=100$ s: step increase in $h_1$ setpoint by +2 cm
$t=300$ s: step increase in $h_2$ setpoint by +2 cm

In [ ]:

Copied!





class SetpointGenerator(Component):
    """Generates step changes in setpoints for h1 and h2."""

    io = IO(outputs=["sp_h1", "sp_h2"])

    def __init__(
        self,
        h1_base: float = H1_0,
        h2_base: float = H2_0,
        h1_step: float = 2.0,
        h2_step: float = 2.0,
        h1_step_time: int = 100,
        h2_step_time: int = 300,
        n_steps: int = N_STEPS,
        **kwargs: _t.Unpack[ComponentArgsDict],
    ) -> None:
        super().__init__(**kwargs)
        self._h1_base = h1_base
        self._h2_base = h2_base
        self._h1_step = h1_step
        self._h2_step = h2_step
        self._h1_step_time = h1_step_time
        self._h2_step_time = h2_step_time
        self._n_steps = n_steps
        self._iteration = 0

    async def step(self) -> None:
        t = self._iteration

        self.sp_h1 = self._h1_base + (self._h1_step if t >= self._h1_step_time else 0.0)
        self.sp_h2 = self._h2_base + (self._h2_step if t >= self._h2_step_time else 0.0)

        self._iteration += 1
        if self._iteration >= self._n_steps:
            await self.io.close()
class SetpointGenerator(Component):
    """Generates step changes in setpoints for h1 and h2."""

    io = IO(outputs=["sp_h1", "sp_h2"])

    def __init__(
        self,
        h1_base: float = H1_0,
        h2_base: float = H2_0,
        h1_step: float = 2.0,
        h2_step: float = 2.0,
        h1_step_time: int = 100,
        h2_step_time: int = 300,
        n_steps: int = N_STEPS,
        **kwargs: _t.Unpack[ComponentArgsDict],
    ) -> None:
        super().__init__(**kwargs)
        self._h1_base = h1_base
        self._h2_base = h2_base
        self._h1_step = h1_step
        self._h2_step = h2_step
        self._h1_step_time = h1_step_time
        self._h2_step_time = h2_step_time
        self._n_steps = n_steps
        self._iteration = 0

    async def step(self) -> None:
        t = self._iteration

        self.sp_h1 = self._h1_base + (self._h1_step if t >= self._h1_step_time else 0.0)
        self.sp_h2 = self._h2_base + (self._h2_step if t >= self._h2_step_time else 0.0)

        self._iteration += 1
        if self._iteration >= self._n_steps:
            await self.io.close()

In [ ]:

Copied!





def build_process(
    gamma1: float,
    gamma2: float,
    kp: tuple[float, float],
    ki: tuple[float, float],
    csv_path: str,
    h1_init: float,
    h2_init: float,
    h3_init: float,
    h4_init: float,
) -> LocalProcess:
    """Factory for a closed-loop quadruple-tank process.

    All connector wiring is identical across scenarios; only the valve ratios,
    controller gains, and initial conditions differ.
    """

    def _connect(src: str, tgt: str) -> AsyncioConnector:
        return AsyncioConnector(spec=ConnectorSpec(source=src, target=tgt))

    comps = [
        SetpointGenerator(name="setpoints"),
        PIController(
            name="controller", kp=kp, ki=ki, initial_values={"h1": [h1_init], "h2": [h2_init]}
        ),
        LowerTank(
            name="tank1",
            tank_area=A1,
            outlet_area=a1,
            upper_outlet_area=a3,
            pump_gain=k1,
            gamma=gamma1,
            initial_values={"prev_h": [h1_init]},
        ),
        LowerTank(
            name="tank2",
            tank_area=A2,
            outlet_area=a2,
            upper_outlet_area=a4,
            pump_gain=k2,
            gamma=gamma2,
            initial_values={"prev_h": [h2_init]},
        ),
        UpperTank(
            name="tank3",
            tank_area=A3,
            outlet_area=a3,
            pump_gain=k2,
            gamma_complement=1.0 - gamma2,
            initial_values={"prev_h": [h3_init]},
        ),
        UpperTank(
            name="tank4",
            tank_area=A4,
            outlet_area=a4,
            pump_gain=k1,
            gamma_complement=1.0 - gamma1,
            initial_values={"prev_h": [h4_init]},
        ),
        FileWriter(
            name="save",
            path=csv_path,
            field_names=["h1", "h2", "h3", "h4", "v1", "v2", "sp_h1", "sp_h2"],
        ),
    ]

    conns = [
        # Setpoints → Controller
        _connect("setpoints.sp_h1", "controller.sp_h1"),
        _connect("setpoints.sp_h2", "controller.sp_h2"),
        # Controller → Pumps → Tanks (direct and cross-coupled flow)
        _connect("controller.v1", "tank1.pump_voltage"),  # Pump 1 → Tank 1 (lower left)
        _connect("controller.v2", "tank2.pump_voltage"),  # Pump 2 → Tank 2 (lower right)
        _connect(
            "controller.v2", "tank3.pump_voltage"
        ),  # Pump 2 → Tank 3 (upper left, cross-coupled)
        _connect(
            "controller.v1", "tank4.pump_voltage"
        ),  # Pump 1 → Tank 4 (upper right, cross-coupled)
        # Upper tanks drain into lower tanks
        _connect("tank3.h", "tank1.upper_h"),
        _connect("tank4.h", "tank2.upper_h"),
        # Tank feedback → Controller
        _connect("tank1.h", "controller.h1"),
        _connect("tank2.h", "controller.h2"),
        # Self-feedback for Euler state (h at step t becomes prev_h at step t+1)
        _connect("tank1.h", "tank1.prev_h"),
        _connect("tank2.h", "tank2.prev_h"),
        _connect("tank3.h", "tank3.prev_h"),
        _connect("tank4.h", "tank4.prev_h"),
        # Save all signals
        _connect("tank1.h", "save.h1"),
        _connect("tank2.h", "save.h2"),
        _connect("tank3.h", "save.h3"),
        _connect("tank4.h", "save.h4"),
        _connect("controller.v1", "save.v1"),
        _connect("controller.v2", "save.v2"),
        _connect("setpoints.sp_h1", "save.sp_h1"),
        _connect("setpoints.sp_h2", "save.sp_h2"),
    ]

    return LocalProcess(components=comps, connectors=conns)
def build_process(
    gamma1: float,
    gamma2: float,
    kp: tuple[float, float],
    ki: tuple[float, float],
    csv_path: str,
    h1_init: float,
    h2_init: float,
    h3_init: float,
    h4_init: float,
) -> LocalProcess:
    """Factory for a closed-loop quadruple-tank process.

    All connector wiring is identical across scenarios; only the valve ratios,
    controller gains, and initial conditions differ.
    """

    def _connect(src: str, tgt: str) -> AsyncioConnector:
        return AsyncioConnector(spec=ConnectorSpec(source=src, target=tgt))

    comps = [
        SetpointGenerator(name="setpoints"),
        PIController(
            name="controller", kp=kp, ki=ki, initial_values={"h1": [h1_init], "h2": [h2_init]}
        ),
        LowerTank(
            name="tank1",
            tank_area=A1,
            outlet_area=a1,
            upper_outlet_area=a3,
            pump_gain=k1,
            gamma=gamma1,
            initial_values={"prev_h": [h1_init]},
        ),
        LowerTank(
            name="tank2",
            tank_area=A2,
            outlet_area=a2,
            upper_outlet_area=a4,
            pump_gain=k2,
            gamma=gamma2,
            initial_values={"prev_h": [h2_init]},
        ),
        UpperTank(
            name="tank3",
            tank_area=A3,
            outlet_area=a3,
            pump_gain=k2,
            gamma_complement=1.0 - gamma2,
            initial_values={"prev_h": [h3_init]},
        ),
        UpperTank(
            name="tank4",
            tank_area=A4,
            outlet_area=a4,
            pump_gain=k1,
            gamma_complement=1.0 - gamma1,
            initial_values={"prev_h": [h4_init]},
        ),
        FileWriter(
            name="save",
            path=csv_path,
            field_names=["h1", "h2", "h3", "h4", "v1", "v2", "sp_h1", "sp_h2"],
        ),
    ]

    conns = [
        # Setpoints → Controller
        _connect("setpoints.sp_h1", "controller.sp_h1"),
        _connect("setpoints.sp_h2", "controller.sp_h2"),
        # Controller → Pumps → Tanks (direct and cross-coupled flow)
        _connect("controller.v1", "tank1.pump_voltage"),  # Pump 1 → Tank 1 (lower left)
        _connect("controller.v2", "tank2.pump_voltage"),  # Pump 2 → Tank 2 (lower right)
        _connect(
            "controller.v2", "tank3.pump_voltage"
        ),  # Pump 2 → Tank 3 (upper left, cross-coupled)
        _connect(
            "controller.v1", "tank4.pump_voltage"
        ),  # Pump 1 → Tank 4 (upper right, cross-coupled)
        # Upper tanks drain into lower tanks
        _connect("tank3.h", "tank1.upper_h"),
        _connect("tank4.h", "tank2.upper_h"),
        # Tank feedback → Controller
        _connect("tank1.h", "controller.h1"),
        _connect("tank2.h", "controller.h2"),
        # Self-feedback for Euler state (h at step t becomes prev_h at step t+1)
        _connect("tank1.h", "tank1.prev_h"),
        _connect("tank2.h", "tank2.prev_h"),
        _connect("tank3.h", "tank3.prev_h"),
        _connect("tank4.h", "tank4.prev_h"),
        # Save all signals
        _connect("tank1.h", "save.h1"),
        _connect("tank2.h", "save.h2"),
        _connect("tank3.h", "save.h3"),
        _connect("tank4.h", "save.h4"),
        _connect("controller.v1", "save.v1"),
        _connect("controller.v2", "save.v2"),
        _connect("setpoints.sp_h1", "save.sp_h1"),
        _connect("setpoints.sp_h2", "save.sp_h2"),
    ]

    return LocalProcess(components=comps, connectors=conns)

Assemble the process¶

Each tank is a separate component, so the physical coupling between them is explicit in the connector wiring.

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.

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.

In [ ]:

Copied!





# Minimum-phase simulation using nominal Johansson (2000) parameter values
process = build_process(
    gamma1=GAMMA1_MP,
    gamma2=GAMMA2_MP,
    kp=(3.0, 2.7),
    ki=(0.1, 0.068),
    csv_path="quadruple_tank_mp.csv",
    h1_init=H1_0,
    h2_init=H2_0,
    h3_init=H3_0,
    h4_init=H4_0,
)
print(f"{len(process.components)} components, {len(process.connectors)} connectors")
# Minimum-phase simulation using nominal Johansson (2000) parameter values
process = build_process(
    gamma1=GAMMA1_MP,
    gamma2=GAMMA2_MP,
    kp=(3.0, 2.7),
    ki=(0.1, 0.068),
    csv_path="quadruple_tank_mp.csv",
    h1_init=H1_0,
    h2_init=H2_0,
    h3_init=H3_0,
    h4_init=H4_0,
)
print(f"{len(process.components)} components, {len(process.connectors)} connectors")

In [ ]:

Copied!

from IPython.display import Image
from plugboard.diagram import MermaidDiagram

Image(url=MermaidDiagram.from_process(process).url)
from IPython.display import Image
from plugboard.diagram import MermaidDiagram

Image(url=MermaidDiagram.from_process(process).url)

Run the simulation¶

In [ ]:

Copied!





async with process:
    await process.run()

df_mp = pd.read_csv("quadruple_tank_mp.csv")
df_mp["time"] = df_mp.index * DT
df_mp.head()
async with process:
    await process.run()

df_mp = pd.read_csv("quadruple_tank_mp.csv")
df_mp["time"] = df_mp.index * DT
df_mp.head()

Plot tank levels over time¶

The top subplot shows the controlled outputs ($h_1$, $h_2$) and their setpoints. The bottom subplot shows the upper tank levels ($h_3$, $h_4$).

In [ ]:

Copied!





fig = make_subplots(
    rows=2,
    cols=1,
    shared_xaxes=True,
    subplot_titles=("Lower tank levels (controlled)", "Upper tank levels"),
    vertical_spacing=0.1,
)

# Lower tank levels and setpoints
fig.add_trace(
    go.Scatter(x=df_mp["time"], y=df_mp["h1"], name="h₁", line=dict(color="blue")), row=1, col=1
)
fig.add_trace(
    go.Scatter(x=df_mp["time"], y=df_mp["h2"], name="h₂", line=dict(color="red")), row=1, col=1
)
fig.add_trace(
    go.Scatter(
        x=df_mp["time"], y=df_mp["sp_h1"], name="SP h₁", line=dict(color="blue", dash="dash")
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=df_mp["time"], y=df_mp["sp_h2"], name="SP h₂", line=dict(color="red", dash="dash")
    ),
    row=1,
    col=1,
)

# Upper tank levels
fig.add_trace(
    go.Scatter(x=df_mp["time"], y=df_mp["h3"], name="h₃", line=dict(color="green")), row=2, col=1
)
fig.add_trace(
    go.Scatter(x=df_mp["time"], y=df_mp["h4"], name="h₄", line=dict(color="orange")), row=2, col=1
)

fig.update_layout(height=600, title_text="Quadruple Tank Process — Minimum Phase (γ₁=0.7, γ₂=0.6)")
fig.update_xaxes(title_text="Time (s)", row=2, col=1)
fig.update_yaxes(title_text="Level (cm)", row=1, col=1)
fig.update_yaxes(title_text="Level (cm)", row=2, col=1)
fig.show()
fig = make_subplots(
    rows=2,
    cols=1,
    shared_xaxes=True,
    subplot_titles=("Lower tank levels (controlled)", "Upper tank levels"),
    vertical_spacing=0.1,
)

# Lower tank levels and setpoints
fig.add_trace(
    go.Scatter(x=df_mp["time"], y=df_mp["h1"], name="h₁", line=dict(color="blue")), row=1, col=1
)
fig.add_trace(
    go.Scatter(x=df_mp["time"], y=df_mp["h2"], name="h₂", line=dict(color="red")), row=1, col=1
)
fig.add_trace(
    go.Scatter(
        x=df_mp["time"], y=df_mp["sp_h1"], name="SP h₁", line=dict(color="blue", dash="dash")
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=df_mp["time"], y=df_mp["sp_h2"], name="SP h₂", line=dict(color="red", dash="dash")
    ),
    row=1,
    col=1,
)

# Upper tank levels
fig.add_trace(
    go.Scatter(x=df_mp["time"], y=df_mp["h3"], name="h₃", line=dict(color="green")), row=2, col=1
)
fig.add_trace(
    go.Scatter(x=df_mp["time"], y=df_mp["h4"], name="h₄", line=dict(color="orange")), row=2, col=1
)

fig.update_layout(height=600, title_text="Quadruple Tank Process — Minimum Phase (γ₁=0.7, γ₂=0.6)")
fig.update_xaxes(title_text="Time (s)", row=2, col=1)
fig.update_yaxes(title_text="Level (cm)", row=1, col=1)
fig.update_yaxes(title_text="Level (cm)", row=2, col=1)
fig.show()

Plot control inputs¶

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).

In [ ]:

Copied!





fig = go.Figure()
fig.add_trace(go.Scatter(x=df_mp["time"], y=df_mp["v1"], name="v₁", line=dict(color="blue")))
fig.add_trace(go.Scatter(x=df_mp["time"], y=df_mp["v2"], name="v₂", line=dict(color="red")))
fig.add_hline(y=0.0, line_dash="dash", line_color="gray", annotation_text="Min (0V)")
fig.add_hline(y=10.0, line_dash="dash", line_color="gray", annotation_text="Max (10V)")
fig.update_layout(
    title="Pump Voltages — Minimum Phase",
    xaxis_title="Time (s)",
    yaxis_title="Voltage (V)",
    height=400,
)
fig.show()
fig = go.Figure()
fig.add_trace(go.Scatter(x=df_mp["time"], y=df_mp["v1"], name="v₁", line=dict(color="blue")))
fig.add_trace(go.Scatter(x=df_mp["time"], y=df_mp["v2"], name="v₂", line=dict(color="red")))
fig.add_hline(y=0.0, line_dash="dash", line_color="gray", annotation_text="Min (0V)")
fig.add_hline(y=10.0, line_dash="dash", line_color="gray", annotation_text="Max (10V)")
fig.update_layout(
    title="Pump Voltages — Minimum Phase",
    xaxis_title="Time (s)",
    yaxis_title="Voltage (V)",
    height=400,
)
fig.show()

Setpoint tracking performance¶

We compute the tracking error for each loop and report key performance metrics.

In [ ]:

Copied!





df_mp["e1"] = df_mp["sp_h1"] - df_mp["h1"]
df_mp["e2"] = df_mp["sp_h2"] - df_mp["h2"]

fig = go.Figure()
fig.add_trace(go.Scatter(x=df_mp["time"], y=df_mp["e1"], name="Error h₁", line=dict(color="blue")))
fig.add_trace(go.Scatter(x=df_mp["time"], y=df_mp["e2"], name="Error h₂", line=dict(color="red")))
fig.add_hline(y=0.0, line_dash="dash", line_color="gray")
fig.update_layout(
    title="Tracking Error — Minimum Phase",
    xaxis_title="Time (s)",
    yaxis_title="Error (cm)",
    height=400,
)
fig.show()
df_mp["e1"] = df_mp["sp_h1"] - df_mp["h1"]
df_mp["e2"] = df_mp["sp_h2"] - df_mp["h2"]

fig = go.Figure()
fig.add_trace(go.Scatter(x=df_mp["time"], y=df_mp["e1"], name="Error h₁", line=dict(color="blue")))
fig.add_trace(go.Scatter(x=df_mp["time"], y=df_mp["e2"], name="Error h₂", line=dict(color="red")))
fig.add_hline(y=0.0, line_dash="dash", line_color="gray")
fig.update_layout(
    title="Tracking Error — Minimum Phase",
    xaxis_title="Time (s)",
    yaxis_title="Error (cm)",
    height=400,
)
fig.show()

In [ ]:

Copied!





def compute_metrics(
    df: pd.DataFrame, sp_col: str, val_col: str, step_time: float
) -> dict[str, float]:
    """Compute step response metrics for a single control loop.

    Parameters
    ----------
    step_time:
        Time (in seconds, matching ``df["time"]``) at which the step change occurs.
    """
    mask = df["time"] >= step_time
    subset = df.loc[mask].copy()

    if subset.empty:
        return {
            "rise_time": float("nan"),
            "settling_time_from_step": float("nan"),
            "overshoot_pct": float("nan"),
            "ss_error": float("nan"),
        }

    try:
        sp_final = float(subset[sp_col].iloc[-1])
    except (KeyError, IndexError, ValueError):
        return {
            "rise_time": float("nan"),
            "settling_time_from_step": float("nan"),
            "overshoot_pct": float("nan"),
            "ss_error": float("nan"),
        }

    try:
        sp_initial = float(df.loc[df["time"] < step_time, sp_col].iloc[-1])
    except (IndexError, KeyError):
        sp_initial = float(subset[sp_col].iloc[0])

    step_size = sp_final - sp_initial

    if abs(step_size) < 1e-6:
        return {
            "rise_time": float("nan"),
            "settling_time_from_step": float("nan"),
            "overshoot_pct": float("nan"),
            "ss_error": float("nan"),
        }

    direction = 1 if step_size > 0 else -1

    val_10 = sp_initial + 0.1 * step_size
    val_90 = sp_initial + 0.9 * step_size

    if direction > 0:
        cond_10 = subset[val_col] >= val_10
        cond_90 = subset[val_col] >= val_90
    else:
        cond_10 = subset[val_col] <= val_10
        cond_90 = subset[val_col] <= val_90

    t_10 = subset.loc[cond_10, "time"].iloc[0] if cond_10.any() else float("nan")
    t_90 = subset.loc[cond_90, "time"].iloc[0] if cond_90.any() else float("nan")

    rise_time = t_90 - t_10 if not (math.isnan(t_10) or math.isnan(t_90)) else float("nan")

    # Settling time: last time error exceeds 2% of step size, measured from the step change
    error = abs(subset[val_col] - sp_final)
    threshold = 0.02 * abs(step_size)
    settled = error <= threshold
    if settled.all():
        settling_time = 0.0
    elif settled.any():
        # last time signal is outside tolerance
        last_out_of_bound = subset.loc[~settled, "time"].iloc[-1]
        settling_time = last_out_of_bound - step_time
    else:
        settling_time = float("nan")

    # Overshoot: consider direction of step
    peak = subset[val_col].max() if direction > 0 else subset[val_col].min()
    if direction > 0:
        overshoot_pct = max(0.0, (peak - sp_final) / abs(step_size) * 100.0)
    else:
        overshoot_pct = max(0.0, (sp_final - peak) / abs(step_size) * 100.0)

    # Steady-state error (average over last 10% of data)
    tail = subset.tail(max(1, len(subset) // 10))
    ss_error = abs(tail[sp_col].mean() - tail[val_col].mean())

    return {
        "rise_time": round(rise_time, 1) if not math.isnan(rise_time) else float("nan"),
        "settling_time_from_step": round(settling_time, 1)
        if not math.isnan(settling_time)
        else float("nan"),
        "overshoot_pct": round(overshoot_pct, 1),
        "ss_error": round(ss_error, 3),
    }


metrics_h1 = compute_metrics(df_mp, "sp_h1", "h1", 100.0)
metrics_h2 = compute_metrics(df_mp, "sp_h2", "h2", 300.0)

pd.DataFrame({"h₁ loop": metrics_h1, "h₂ loop": metrics_h2}).rename(
    index={
        "rise_time": "Rise time (s)",
        "settling_time_from_step": "Settling time from step (s)",
        "overshoot_pct": "Overshoot (%)",
        "ss_error": "Steady-state error (cm)",
    },
)
def compute_metrics(
    df: pd.DataFrame, sp_col: str, val_col: str, step_time: float
) -> dict[str, float]:
    """Compute step response metrics for a single control loop.

    Parameters
    ----------
    step_time:
        Time (in seconds, matching ``df["time"]``) at which the step change occurs.
    """
    mask = df["time"] >= step_time
    subset = df.loc[mask].copy()

    if subset.empty:
        return {
            "rise_time": float("nan"),
            "settling_time_from_step": float("nan"),
            "overshoot_pct": float("nan"),
            "ss_error": float("nan"),
        }

    try:
        sp_final = float(subset[sp_col].iloc[-1])
    except (KeyError, IndexError, ValueError):
        return {
            "rise_time": float("nan"),
            "settling_time_from_step": float("nan"),
            "overshoot_pct": float("nan"),
            "ss_error": float("nan"),
        }

    try:
        sp_initial = float(df.loc[df["time"] < step_time, sp_col].iloc[-1])
    except (IndexError, KeyError):
        sp_initial = float(subset[sp_col].iloc[0])

    step_size = sp_final - sp_initial

    if abs(step_size) < 1e-6:
        return {
            "rise_time": float("nan"),
            "settling_time_from_step": float("nan"),
            "overshoot_pct": float("nan"),
            "ss_error": float("nan"),
        }

    direction = 1 if step_size > 0 else -1

    val_10 = sp_initial + 0.1 * step_size
    val_90 = sp_initial + 0.9 * step_size

    if direction > 0:
        cond_10 = subset[val_col] >= val_10
        cond_90 = subset[val_col] >= val_90
    else:
        cond_10 = subset[val_col] <= val_10
        cond_90 = subset[val_col] <= val_90

    t_10 = subset.loc[cond_10, "time"].iloc[0] if cond_10.any() else float("nan")
    t_90 = subset.loc[cond_90, "time"].iloc[0] if cond_90.any() else float("nan")

    rise_time = t_90 - t_10 if not (math.isnan(t_10) or math.isnan(t_90)) else float("nan")

    # Settling time: last time error exceeds 2% of step size, measured from the step change
    error = abs(subset[val_col] - sp_final)
    threshold = 0.02 * abs(step_size)
    settled = error <= threshold
    if settled.all():
        settling_time = 0.0
    elif settled.any():
        # last time signal is outside tolerance
        last_out_of_bound = subset.loc[~settled, "time"].iloc[-1]
        settling_time = last_out_of_bound - step_time
    else:
        settling_time = float("nan")

    # Overshoot: consider direction of step
    peak = subset[val_col].max() if direction > 0 else subset[val_col].min()
    if direction > 0:
        overshoot_pct = max(0.0, (peak - sp_final) / abs(step_size) * 100.0)
    else:
        overshoot_pct = max(0.0, (sp_final - peak) / abs(step_size) * 100.0)

    # Steady-state error (average over last 10% of data)
    tail = subset.tail(max(1, len(subset) // 10))
    ss_error = abs(tail[sp_col].mean() - tail[val_col].mean())

    return {
        "rise_time": round(rise_time, 1) if not math.isnan(rise_time) else float("nan"),
        "settling_time_from_step": round(settling_time, 1)
        if not math.isnan(settling_time)
        else float("nan"),
        "overshoot_pct": round(overshoot_pct, 1),
        "ss_error": round(ss_error, 3),
    }


metrics_h1 = compute_metrics(df_mp, "sp_h1", "h1", 100.0)
metrics_h2 = compute_metrics(df_mp, "sp_h2", "h2", 300.0)

pd.DataFrame({"h₁ loop": metrics_h1, "h₂ loop": metrics_h2}).rename(
    index={
        "rise_time": "Rise time (s)",
        "settling_time_from_step": "Settling time from step (s)",
        "overshoot_pct": "Overshoot (%)",
        "ss_error": "Steady-state error (cm)",
    },
)

Minimum-phase vs non-minimum-phase¶

The key insight of the quadruple tank process is that valve parameters $\gamma_1$ and $\gamma_2$ determine the system's phase characteristics:

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.
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.

Let's re-run the simulation with non-minimum-phase parameters and compare.

In [ ]:

Copied!





# Non-minimum-phase simulation — uses correct NMP steady-state initial conditions
# and more conservative PI gains (NMP systems require reduced bandwidth)
process_nmp = build_process(
    gamma1=GAMMA1_NMP,
    gamma2=GAMMA2_NMP,
    kp=(0.5, 0.5),
    ki=(0.005, 0.004),
    csv_path="quadruple_tank_nmp.csv",
    h1_init=H1_NMP,
    h2_init=H2_NMP,
    h3_init=H3_NMP,
    h4_init=H4_NMP,
)

async with process_nmp:
    await process_nmp.run()

df_nmp = pd.read_csv("quadruple_tank_nmp.csv")
df_nmp["time"] = df_nmp.index * DT
# Non-minimum-phase simulation — uses correct NMP steady-state initial conditions
# and more conservative PI gains (NMP systems require reduced bandwidth)
process_nmp = build_process(
    gamma1=GAMMA1_NMP,
    gamma2=GAMMA2_NMP,
    kp=(0.5, 0.5),
    ki=(0.005, 0.004),
    csv_path="quadruple_tank_nmp.csv",
    h1_init=H1_NMP,
    h2_init=H2_NMP,
    h3_init=H3_NMP,
    h4_init=H4_NMP,
)

async with process_nmp:
    await process_nmp.run()

df_nmp = pd.read_csv("quadruple_tank_nmp.csv")
df_nmp["time"] = df_nmp.index * DT

In [ ]:

Copied!





fig = make_subplots(
    rows=2,
    cols=2,
    shared_xaxes=True,
    subplot_titles=(
        "Minimum phase (γ₁=0.7, γ₂=0.6)",
        "Non-minimum phase (γ₁=0.43, γ₂=0.34)",
        "Upper tanks (MP)",
        "Upper tanks (NMP)",
    ),
    vertical_spacing=0.12,
    horizontal_spacing=0.08,
)

# Minimum phase — lower tanks
fig.add_trace(
    go.Scatter(
        x=df_mp["time"], y=df_mp["h1"], name="h₁ (MP)", line=dict(color="blue"), legendgroup="mp"
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=df_mp["time"], y=df_mp["h2"], name="h₂ (MP)", line=dict(color="red"), legendgroup="mp"
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=df_mp["time"],
        y=df_mp["sp_h1"],
        name="SP h₁",
        line=dict(color="blue", dash="dash"),
        legendgroup="mp",
        showlegend=False,
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=df_mp["time"],
        y=df_mp["sp_h2"],
        name="SP h₂",
        line=dict(color="red", dash="dash"),
        legendgroup="mp",
        showlegend=False,
    ),
    row=1,
    col=1,
)

# Non-minimum phase — lower tanks
fig.add_trace(
    go.Scatter(
        x=df_nmp["time"],
        y=df_nmp["h1"],
        name="h₁ (NMP)",
        line=dict(color="blue", dash="dot"),
        legendgroup="nmp",
    ),
    row=1,
    col=2,
)
fig.add_trace(
    go.Scatter(
        x=df_nmp["time"],
        y=df_nmp["h2"],
        name="h₂ (NMP)",
        line=dict(color="red", dash="dot"),
        legendgroup="nmp",
    ),
    row=1,
    col=2,
)
fig.add_trace(
    go.Scatter(
        x=df_nmp["time"],
        y=df_nmp["sp_h1"],
        name="SP h₁",
        line=dict(color="blue", dash="dash"),
        legendgroup="nmp",
        showlegend=False,
    ),
    row=1,
    col=2,
)
fig.add_trace(
    go.Scatter(
        x=df_nmp["time"],
        y=df_nmp["sp_h2"],
        name="SP h₂",
        line=dict(color="red", dash="dash"),
        legendgroup="nmp",
        showlegend=False,
    ),
    row=1,
    col=2,
)

# Minimum phase — upper tanks
fig.add_trace(
    go.Scatter(
        x=df_mp["time"], y=df_mp["h3"], name="h₃ (MP)", line=dict(color="green"), legendgroup="mp"
    ),
    row=2,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=df_mp["time"], y=df_mp["h4"], name="h₄ (MP)", line=dict(color="orange"), legendgroup="mp"
    ),
    row=2,
    col=1,
)

# Non-minimum phase — upper tanks
fig.add_trace(
    go.Scatter(
        x=df_nmp["time"],
        y=df_nmp["h3"],
        name="h₃ (NMP)",
        line=dict(color="green", dash="dot"),
        legendgroup="nmp",
    ),
    row=2,
    col=2,
)
fig.add_trace(
    go.Scatter(
        x=df_nmp["time"],
        y=df_nmp["h4"],
        name="h₄ (NMP)",
        line=dict(color="orange", dash="dot"),
        legendgroup="nmp",
    ),
    row=2,
    col=2,
)

fig.update_layout(height=700, title_text="Minimum Phase vs Non-Minimum Phase Comparison")
fig.update_xaxes(title_text="Time (s)", row=2, col=1)
fig.update_xaxes(title_text="Time (s)", row=2, col=2)
for row in [1, 2]:
    for col in [1, 2]:
        fig.update_yaxes(title_text="Level (cm)", row=row, col=col)
fig.show()
fig = make_subplots(
    rows=2,
    cols=2,
    shared_xaxes=True,
    subplot_titles=(
        "Minimum phase (γ₁=0.7, γ₂=0.6)",
        "Non-minimum phase (γ₁=0.43, γ₂=0.34)",
        "Upper tanks (MP)",
        "Upper tanks (NMP)",
    ),
    vertical_spacing=0.12,
    horizontal_spacing=0.08,
)

# Minimum phase — lower tanks
fig.add_trace(
    go.Scatter(
        x=df_mp["time"], y=df_mp["h1"], name="h₁ (MP)", line=dict(color="blue"), legendgroup="mp"
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=df_mp["time"], y=df_mp["h2"], name="h₂ (MP)", line=dict(color="red"), legendgroup="mp"
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=df_mp["time"],
        y=df_mp["sp_h1"],
        name="SP h₁",
        line=dict(color="blue", dash="dash"),
        legendgroup="mp",
        showlegend=False,
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=df_mp["time"],
        y=df_mp["sp_h2"],
        name="SP h₂",
        line=dict(color="red", dash="dash"),
        legendgroup="mp",
        showlegend=False,
    ),
    row=1,
    col=1,
)

# Non-minimum phase — lower tanks
fig.add_trace(
    go.Scatter(
        x=df_nmp["time"],
        y=df_nmp["h1"],
        name="h₁ (NMP)",
        line=dict(color="blue", dash="dot"),
        legendgroup="nmp",
    ),
    row=1,
    col=2,
)
fig.add_trace(
    go.Scatter(
        x=df_nmp["time"],
        y=df_nmp["h2"],
        name="h₂ (NMP)",
        line=dict(color="red", dash="dot"),
        legendgroup="nmp",
    ),
    row=1,
    col=2,
)
fig.add_trace(
    go.Scatter(
        x=df_nmp["time"],
        y=df_nmp["sp_h1"],
        name="SP h₁",
        line=dict(color="blue", dash="dash"),
        legendgroup="nmp",
        showlegend=False,
    ),
    row=1,
    col=2,
)
fig.add_trace(
    go.Scatter(
        x=df_nmp["time"],
        y=df_nmp["sp_h2"],
        name="SP h₂",
        line=dict(color="red", dash="dash"),
        legendgroup="nmp",
        showlegend=False,
    ),
    row=1,
    col=2,
)

# Minimum phase — upper tanks
fig.add_trace(
    go.Scatter(
        x=df_mp["time"], y=df_mp["h3"], name="h₃ (MP)", line=dict(color="green"), legendgroup="mp"
    ),
    row=2,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=df_mp["time"], y=df_mp["h4"], name="h₄ (MP)", line=dict(color="orange"), legendgroup="mp"
    ),
    row=2,
    col=1,
)

# Non-minimum phase — upper tanks
fig.add_trace(
    go.Scatter(
        x=df_nmp["time"],
        y=df_nmp["h3"],
        name="h₃ (NMP)",
        line=dict(color="green", dash="dot"),
        legendgroup="nmp",
    ),
    row=2,
    col=2,
)
fig.add_trace(
    go.Scatter(
        x=df_nmp["time"],
        y=df_nmp["h4"],
        name="h₄ (NMP)",
        line=dict(color="orange", dash="dot"),
        legendgroup="nmp",
    ),
    row=2,
    col=2,
)

fig.update_layout(height=700, title_text="Minimum Phase vs Non-Minimum Phase Comparison")
fig.update_xaxes(title_text="Time (s)", row=2, col=1)
fig.update_xaxes(title_text="Time (s)", row=2, col=2)
for row in [1, 2]:
    for col in [1, 2]:
        fig.update_yaxes(title_text="Level (cm)", row=row, col=col)
fig.show()

In [ ]:

Copied!





# Zoom in on the NMP h1 response to reveal the inverse‑response (wrong‑direction) transient.
# For a non-minimum-phase system the step in sp_h2 at t=300 s causes h1 to dip *down*
# before recovering — a hallmark characteristic of cross-coupled NMP dynamics.
ZOOM_START, ZOOM_END = 295, 360

mask = (df_nmp["time"] >= ZOOM_START) & (df_nmp["time"] <= ZOOM_END)
zoom = df_nmp.loc[mask]

fig_zoom = go.Figure()
fig_zoom.add_trace(
    go.Scatter(
        x=zoom["time"],
        y=zoom["h1"],
        name="h₁ (NMP)",
        line=dict(color="blue"),
    )
)
fig_zoom.add_trace(
    go.Scatter(
        x=zoom["time"],
        y=zoom["sp_h1"],
        name="SP h₁",
        line=dict(color="blue", dash="dash"),
    )
)
# Mark the step change
fig_zoom.add_vline(
    x=300,
    line_dash="dot",
    line_color="gray",
    annotation_text="sp_h2 step",
    annotation_position="top right",
)
# Annotate the inverse-response dip
inv_idx = zoom["h1"].idxmin()
fig_zoom.add_annotation(
    x=zoom.loc[inv_idx, "time"],
    y=zoom.loc[inv_idx, "h1"],
    text="Inverse response<br>(initial decrease)",
    showarrow=True,
    arrowhead=2,
    ax=40,
    ay=-40,
    font=dict(size=12),
)
fig_zoom.update_layout(
    title="NMP inverse response — zoom around h₂ step change (t = 300 s)",
    xaxis_title="Time (s)",
    yaxis_title="Tank level (cm)",
    legend=dict(orientation="h", y=-0.2),
    template="plotly_white",
)
fig_zoom.show()
# Zoom in on the NMP h1 response to reveal the inverse‑response (wrong‑direction) transient.
# For a non-minimum-phase system the step in sp_h2 at t=300 s causes h1 to dip *down*
# before recovering — a hallmark characteristic of cross-coupled NMP dynamics.
ZOOM_START, ZOOM_END = 295, 360

mask = (df_nmp["time"] >= ZOOM_START) & (df_nmp["time"] <= ZOOM_END)
zoom = df_nmp.loc[mask]

fig_zoom = go.Figure()
fig_zoom.add_trace(
    go.Scatter(
        x=zoom["time"],
        y=zoom["h1"],
        name="h₁ (NMP)",
        line=dict(color="blue"),
    )
)
fig_zoom.add_trace(
    go.Scatter(
        x=zoom["time"],
        y=zoom["sp_h1"],
        name="SP h₁",
        line=dict(color="blue", dash="dash"),
    )
)
# Mark the step change
fig_zoom.add_vline(
    x=300,
    line_dash="dot",
    line_color="gray",
    annotation_text="sp_h2 step",
    annotation_position="top right",
)
# Annotate the inverse-response dip
inv_idx = zoom["h1"].idxmin()
fig_zoom.add_annotation(
    x=zoom.loc[inv_idx, "time"],
    y=zoom.loc[inv_idx, "h1"],
    text="Inverse response<br>(initial decrease)",
    showarrow=True,
    arrowhead=2,
    ax=40,
    ay=-40,
    font=dict(size=12),
)
fig_zoom.update_layout(
    title="NMP inverse response — zoom around h₂ step change (t = 300 s)",
    xaxis_title="Time (s)",
    yaxis_title="Tank level (cm)",
    legend=dict(orientation="h", y=-0.2),
    template="plotly_white",
)
fig_zoom.show()

In [ ]:

Copied!





# Compare performance metrics
metrics_nmp_h1 = compute_metrics(df_nmp, "sp_h1", "h1", 100.0)
metrics_nmp_h2 = compute_metrics(df_nmp, "sp_h2", "h2", 300.0)

pd.DataFrame(
    {
        "h₁ (MP)": metrics_h1,
        "h₁ (NMP)": metrics_nmp_h1,
        "h₂ (MP)": metrics_h2,
        "h₂ (NMP)": metrics_nmp_h2,
    }
).rename(
    index={
        "rise_time": "Rise time (s)",
        "settling_time_from_step": "Settling time from step (s)",
        "overshoot_pct": "Overshoot (%)",
        "ss_error": "Steady-state error (cm)",
    }
)
# Compare performance metrics
metrics_nmp_h1 = compute_metrics(df_nmp, "sp_h1", "h1", 100.0)
metrics_nmp_h2 = compute_metrics(df_nmp, "sp_h2", "h2", 300.0)

pd.DataFrame(
    {
        "h₁ (MP)": metrics_h1,
        "h₁ (NMP)": metrics_nmp_h1,
        "h₂ (MP)": metrics_h2,
        "h₂ (NMP)": metrics_nmp_h2,
    }
).rename(
    index={
        "rise_time": "Rise time (s)",
        "settling_time_from_step": "Settling time from step (s)",
        "overshoot_pct": "Overshoot (%)",
        "ss_error": "Steady-state error (cm)",
    }
)

Summary¶

This notebook demonstrated how to model and simulate the quadruple tank process using Plugboard:

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.
Controller: A decentralized PI controller with anti-windup clamping regulates the lower tank levels.
Setpoint generator: Step changes in both setpoints showcase the MIMO coupling between loops.
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.

Key takeaways¶

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.
The minimum-phase configuration responds faster because pump flow goes directly to the controlled (lower) tanks.
The non-minimum-phase configuration requires more conservative tuning and responds much more slowly, because the flow must pass through the upper tanks first.
This benchmark illustrates why MIMO systems with right-half-plane zeros pose fundamental performance limitations for feedback control.

References¶

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.