Dynamics & process control¶

The differentiable ODE integrators, dynamic holdups and dynamic flowsheets, the classical PID controller, and the tuning/optimal-control helpers for the time domain.

See the dynamics & process control guide for worked examples.

Dynamics¶

dynamics ¶

Differentiable dynamic (time-domain) process simulation for Fugacio.

The fugacio.sim blocks elsewhere are steady-state; this subpackage adds the time dimension while keeping everything end-to-end differentiable:

integrate: the ODE integration core: a fixed output-grid odeint (explicit and stiff implicit steppers, differentiable in both modes) and an adaptive integrate with a continuous-adjoint custom_vjp;
units: dynamic unit operations carried as ODEs in conserved holdups (buffer/level tanks, mixing tanks, dynamic CSTR, dynamic flash, lumped thermal mass / heat exchanger, gas receiver), reusing the steady-state thermodynamics for the instantaneous constitutive relations;
flowsheet: DynamicFlowsheet, a declarative assembly of dynamic units and controllers into one global ODE that is simulated over time and differentiated through;
optimize: dynamic optimization and estimation (optimal control over input trajectories, time-series parameter estimation), composing the existing optimizers with the integrator.

Modules:

Name	Description
`flowsheet`	Assemble dynamic units and control loops into one differentiable flowsheet ODE.
`integrate`	Differentiable time integration of ordinary differential equations.
`optimize`	Dynamic optimization, estimation, and controller tuning.
`units`	Dynamic unit operations: holdups carried as ODE states.

Classes:

Name	Description
`DynamicFlowsheet`	A declarative dynamic flowsheet: feeds, dynamic units, and control loops.
`DynamicResult`	Trajectory returned by `DynamicFlowsheet.simulate`.
`ODEResult`	Outcome of an adaptive `integrate` call.
`DynamicEstimateResult`	Outcome of `estimate_dynamics`.
`OptimalControlResult`	Outcome of `optimal_control`.
`DynamicCSTR`	A non-isothermal, constant-volume liquid CSTR: the canonical control plant.
`DynamicFlash`	An isothermal-isobaric flash drum with liquid holdup and equilibrium vapor.
`DynamicUnit`	Abstract base class for a dynamic unit operation.
`GasReceiver`	An isothermal gas surge drum: pressure builds and decays with the inventory.
`LevelTank`	An isothermal, well-mixed liquid surge tank with level and composition dynamics.
`MixingTank`	A constant-holdup, perfectly mixed blending tank (composition first-order lag).
`ThermalMass`	A heated/cooled, well-mixed stirred tank: lumped temperature dynamics.

Functions:

Name	Description
`simulate`	Integrate a bare `rhs(t, y, theta)` over `ts` (a thin `odeint` alias).
`odeint`	Integrate `dy/dt = func(t, y, theta)` over the output grid `ts`.
`odeint_final`	Integrate from `t0` to `t1` and return only the final state.
`estimate_dynamics`	Fit parameters `theta` so the simulated trajectory matches `data`.
`optimal_control`	Find the input trajectory minimizing a running-plus-terminal cost.
`tune_pid`	Tune controller gains by minimizing a closed-loop error integral.

DynamicFlowsheet `dataclass` ¶

DynamicFlowsheet(
    feeds: dict[str, FeedSpec] = dict(),
    units: dict[str, _Unit] = dict(),
    loops: list[_Loop] = list(),
    manips: list[_Manip] = list(),
    _order: list[str] = list(),
)

A declarative dynamic flowsheet: feeds, dynamic units, and control loops.

Example::

fs = DynamicFlowsheet()
fs.feed("feed", feed_stream)                       # constant or (t, theta) -> Stream
fs.add(tank, inputs=["feed"])                      # a LevelTank named "tank"
fs.control(pid, measurement=("tank", "level"),     # level controller ...
           setpoint=2.0, actuator=("tank", "flow"))# ... manipulating outlet flow
result = fs.simulate(ts=jnp.linspace(0, 100, 201))
level = result.measurement("tank", "level")

Units are advanced in registration order; reference an upstream outlet as "unit.port" (port index) and a feed by its name.

Methods:

Name	Description
`feed`	Register a feed stream (a fixed `Stream` or `(t, theta) -> Stream`).
`add`	Register a dynamic `unit` with its inlet sources (feeds or `"unit.port"`).
`control`	Close a feedback loop: `pid` drives `actuator` to hold `measurement`.
`manipulate`	Drive a manipulated variable open-loop from a schedule `fn(t, theta)`.
`initial_state`	Build the global initial state (unit holdups + controller states).
`rhs`	Global state derivative `d(state)/dt` (the function handed to the integrator).
`simulate`	Integrate the flowsheet over the output grid `ts` and return a `DynamicResult`.

feed ¶

feed(name: str, spec: FeedSpec) -> DynamicFlowsheet

Register a feed stream (a fixed Stream or (t, theta) -> Stream).

add ¶

add(
    unit: DynamicUnit, *, inputs: Sequence[str] = ()
) -> DynamicFlowsheet

Register a dynamic unit with its inlet sources (feeds or "unit.port").

control ¶

control(
    pid: PID,
    *,
    measurement: MeasurementSpec,
    setpoint: SetpointSpec,
    actuator: tuple[str, str],
    u0: ArrayLike | None = None,
) -> DynamicFlowsheet

Close a feedback loop: pid drives actuator to hold measurement.

manipulate ¶

manipulate(
    unit: str, mv: str, fn: Callable[[Array, Any], Array]
) -> DynamicFlowsheet

Drive a manipulated variable open-loop from a schedule fn(t, theta).

initial_state ¶

initial_state(
    t0: ArrayLike = 0.0, theta: Any = None
) -> dict[str, Any]

Build the global initial state (unit holdups + controller states).

rhs ¶

rhs(
    t: Array, state: dict[str, Any], theta: Any = None
) -> dict[str, Any]

Global state derivative d(state)/dt (the function handed to the integrator).

simulate ¶

simulate(
    ts: Array,
    *,
    y0: dict[str, Any] | None = None,
    theta: Any = None,
    method: str = "rk4",
    substeps: int = 4,
) -> DynamicResult

Integrate the flowsheet over the output grid ts and return a DynamicResult.

Parameters:

Name	Type	Description	Default
`ts`	`Array`	Strictly increasing output times.	required
`y0`	`dict[str, Any] \| None`	Optional global initial state (defaults to each unit's initial state and a bumpless controller start).	`None`
`theta`	`Any`	Differentiable parameter pytree threaded to feeds, setpoints, manipulations and units.	`None`
`method`	`str`	Integration method (see `fugacio.sim.dynamics.FIXED_METHODS`).	`'rk4'`
`substeps`	`int`	Inner steps per output interval.	`4`

The returned trajectories of measurements and controller outputs are recomputed from the state trajectory, so they are differentiable in theta and y0 as well.

DynamicResult ¶

Bases: NamedTuple

Trajectory returned by DynamicFlowsheet.simulate.

Attributes:

Name	Type	Description
`ts`	`Array`	Output times (shape `(T,)`).
`states`	`dict[str, Any]`	Per-unit state trajectory (a dict `name -> state` with a leading time axis on every leaf).
`measurements`	`dict[str, dict[str, Array]]`	Per-unit measurement trajectory (a dict `name -> {key: ...}`).
`controls`	`dict[str, Array]`	Manipulated-variable trajectory (a dict `"unit.mv" -> values`).

Methods:

Name	Description
`measurement`	Trajectory of a named measurement for a unit.

measurement ¶

measurement(unit: str, key: str) -> Array

Trajectory of a named measurement for a unit.

ODEResult ¶

Bases: NamedTuple

Outcome of an adaptive integrate call.

Attributes:

Name	Type	Description
`y`	`Any`	Final state pytree at `t1` (differentiable w.r.t. `y0` and `theta`).
`n_steps`	`Array`	Total attempted steps (accepted + rejected).
`n_accepted`	`Array`	Accepted steps.
`success`	`Array`	Whether the march reached `t1` within `max_steps`.

DynamicEstimateResult ¶

Bases: NamedTuple

Outcome of estimate_dynamics.

Attributes:

Name	Type	Description
`theta`	`Any`	Estimated parameter pytree.
`trajectory`	`Any`	Model trajectory at the estimate (leading time axis).
`cost`	`Array`	Half-sum-of-squares residual at the estimate.
`result`	`OptimizeResult`	The underlying `fugacio.sim.OptimizeResult`.

OptimalControlResult ¶

Bases: NamedTuple

Outcome of optimal_control.

Attributes:

Name	Type	Description
`u`	`Any`	Optimal piecewise-constant control (one value, or vector, per interval).
`cost`	`Array`	Optimal total cost (running integral plus terminal).
`trajectory`	`Any`	State trajectory under the optimal control (leading time axis).
`result`	`OptimizeResult`	The underlying `fugacio.sim.OptimizeResult`.

DynamicCSTR `dataclass` ¶

DynamicCSTR(
    name: str,
    components: tuple[str, ...],
    reactions: Reaction | Sequence[Reaction],
    rate_laws: Any,
    volume: ArrayLike = 1.0,
    q: ArrayLike = 1.0,
    rho_cp: ArrayLike = 4000000.0,
    ua: ArrayLike = 0.0,
    jacket_t: ArrayLike = 298.15,
    p: ArrayLike = 101325.0,
    c0: Array | None = None,
    t0: ArrayLike = 298.15,
)

Bases: DynamicUnit

A non-isothermal, constant-volume liquid CSTR: the canonical control plant.

Works in concentration space: state is the vector of species concentrations C (mol/m^3) and the reactor temperature T. The reactor has volume volume and a constant volumetric throughput q, so the residence time is volume / q. The energy balance carries the heat of reaction (from fugacio.thermo.delta_h_rxn) and a jacket term UA (T_jacket - T) with a constant volumetric heat capacity rho_cp (J/m^3/K). Exothermic operation reproduces the classic ignition/extinction and limit-cycle behaviour, which is why this is the reactor-control benchmark.

Manipulated variables: "jacket_t" (K) and "q" (m^3/s). Measurements: temperature, concentration, heat_release.

DynamicFlash `dataclass` ¶

DynamicFlash(
    name: str,
    components: tuple[str, ...],
    t: ArrayLike = 298.15,
    p: ArrayLike = 101325.0,
    eos: CubicEOS = PR,
    kij: Array | None = None,
    vapor_draw: ArrayLike = 0.0,
    liquid_draw: ArrayLike = 0.0,
    m0: Array | None = None,
)

Bases: DynamicUnit

An isothermal-isobaric flash drum with liquid holdup and equilibrium vapor.

State is the per-component liquid holdup M (mol). The vapour drawn off is in instantaneous phase equilibrium with the well-mixed holdup (its composition is the equilibrium vapour of a flash_pt on the holdup at (T, P)) and leaves at a commanded rate; the liquid product leaves at the holdup composition. This captures the composition response of a separator to feed and draw disturbances while reusing the rigorous EOS equilibrium.

Manipulated variables: "vapor_draw" (mol/s), "liquid_draw" (mol/s). Measurements: holdup, x (liquid composition), y (vapour composition).

DynamicUnit ¶

Bases: ABC

Abstract base class for a dynamic unit operation.

Methods:

Name	Description
`initial_state`	Return the initial state pytree for this unit.
`evaluate`	Return the state derivative, outlet streams and measurements at `state`.
`measured`	State-only measurements (level, T, P, composition) available to controllers.

initial_state `abstractmethod` ¶

initial_state() -> Any

Return the initial state pytree for this unit.

evaluate `abstractmethod` ¶

evaluate(
    state: Any,
    inlets: tuple[Stream, ...],
    controls: Controls | None = None,
    theta: Any = None,
) -> UnitStep

Return the state derivative, outlet streams and measurements at state.

measured `abstractmethod` ¶

measured(state: Any) -> dict[str, Array]

State-only measurements (level, T, P, composition) available to controllers.

These depend on the unit state alone, not on the inlet streams, so a controller can read them before the units are advanced, which keeps the instantaneous control algebra explicit (see DynamicFlowsheet).

GasReceiver `dataclass` ¶

GasReceiver(
    name: str,
    components: tuple[str, ...],
    volume: ArrayLike = 1.0,
    t: ArrayLike = 298.15,
    outlet: str = "valve",
    valve_cv: ArrayLike = 0.001,
    p_down: ArrayLike = 101325.0,
    flow_setpoint: ArrayLike = 0.0,
    n0: Array | None = None,
)

Bases: DynamicUnit

An isothermal gas surge drum: pressure builds and decays with the inventory.

State is the component mole holdup N of gas in a fixed volume; the pressure follows the ideal-gas law P = (sum N) R T / V. Gas leaves either at a commanded flow (outlet="flow", manipulated variable "flow") or through a valve to a downstream pressure p_down (outlet="valve", manipulated variable "opening"), giving the textbook pressure-control plant.

Manipulated variables: "flow" (mol/s) or "opening" (-). Measurements: pressure, holdup, outlet_flow.

LevelTank `dataclass` ¶

LevelTank(
    name: str,
    components: tuple[str, ...],
    area: ArrayLike = 1.0,
    t: ArrayLike = 298.15,
    p: ArrayLike = 101325.0,
    outlet: str = "pump",
    valve_cv: ArrayLike = 1.0,
    flow_setpoint: ArrayLike = 0.0,
    n0: Array | None = None,
)

Bases: DynamicUnit

An isothermal, well-mixed liquid surge tank with level and composition dynamics.

State is the per-component liquid holdup N (mol). The liquid leaves either at a commanded molar flow (outlet="pump", manipulated variable "flow") or through a valve whose discharge follows a Torricelli law (outlet="valve", manipulated variable "opening" in [0, 1]); the outlet always carries the well-mixed holdup composition. The liquid level is inferred from the holdup volume and the tank cross-sectional area using the real liquid density.

Manipulated variables: "flow" (mol/s, pump mode) or "opening" (-, valve mode). Measurements: level, volume, holdup, outlet_flow.

MixingTank `dataclass` ¶

MixingTank(
    name: str,
    components: tuple[str, ...],
    holdup: ArrayLike = 1.0,
    t: ArrayLike = 298.15,
    p: ArrayLike = 101325.0,
    x0: Array | None = None,
)

Bases: DynamicUnit

A constant-holdup, perfectly mixed blending tank (composition first-order lag).

With a fixed molar holdup holdup and equal in/out flow, the outlet composition relaxes toward the (instantaneous) inlet composition with a time constant equal to the residence time holdup / inlet_flow. State is the mole-fraction vector x.

Measurements: residence_time, inlet_flow, x0 ... (per-component fractions are also returned under frac).

ThermalMass `dataclass` ¶

ThermalMass(
    name: str,
    components: tuple[str, ...],
    holdup: ArrayLike = 100.0,
    ua: ArrayLike = 0.0,
    t_ambient: ArrayLike = 298.15,
    p: ArrayLike = 101325.0,
    duty_setpoint: ArrayLike = 0.0,
    t0: ArrayLike = 298.15,
)

Bases: DynamicUnit

A heated/cooled, well-mixed stirred tank: lumped temperature dynamics.

A fixed molar holdup of liquid is heated by a duty Q (manipulated variable "duty", W), exchanges sensible heat with the inlet flow, and loses heat to ambient through ua (W/K). State is the bulk temperature T; the outlet carries the inlet flows at T. Heat capacities come from the ideal-gas correlations of fugacio.thermo.

Manipulated variable: "duty" (W). Measurements: temperature, duty.

simulate ¶

simulate(
    rhs: Callable[[Array, Any, Any], Any],
    y0: Any,
    ts: Array,
    theta: Any = None,
    *,
    method: str = "rk4",
    substeps: int = 4,
) -> Any

Integrate a bare rhs(t, y, theta) over ts (a thin odeint alias).

For ad-hoc dynamic models that are not expressed as a DynamicFlowsheet. Returns the state trajectory pytree (leading time axis), differentiable in y0 and theta.

odeint ¶

odeint(
    func: RHS,
    y0: Any,
    ts: Array,
    theta: Any = None,
    *,
    method: str = "rk4",
    substeps: int = 1,
) -> Any

Integrate dy/dt = func(t, y, theta) over the output grid ts.

The state is advanced from ts[0] to ts[-1], taking substeps uniform inner steps of the chosen method between successive output points, and the state is recorded at every entry of ts. Because the step pattern is static the whole integration is an ordinary jax.lax.scan, hence differentiable in forward and reverse mode with respect to y0 and theta, no custom rule needed.

Parameters:

Name	Type	Description	Default
`func`	`RHS`	Right-hand side `func(t, y, theta) -> dy` (`dy` matches `y`'s pytree structure).	required
`y0`	`Any`	Initial state pytree at `ts[0]`.	required
`ts`	`Array`	1-D array of strictly increasing output times (length >= 2).	required
`theta`	`Any`	Optional differentiable parameter pytree forwarded to `func`.	`None`
`method`	`str`	One of `FIXED_METHODS`: `"euler"`, `"rk4"` (default), `"dopri5"`, or the stiff `"implicit_euler"` / `"trapezoidal"`.	`'rk4'`
`substeps`	`int`	Number of inner integration steps per output interval (>= 1); raise it to cut discretisation error without densifying `ts`.	`1`

Returns:

Type	Description
`Any`	The trajectory as a pytree of the same structure as `y0` with a leading
`Any`	time axis of length `len(ts)` on each leaf (so `out[k]` is the state at
`Any`	`ts[k]`). Differentiable with respect to `y0` and `theta`.

odeint_final ¶

odeint_final(
    func: RHS,
    y0: Any,
    t0: ArrayLike,
    t1: ArrayLike,
    theta: Any = None,
    *,
    method: str = "rk4",
    steps: int = 100,
) -> Any

Integrate from t0 to t1 and return only the final state.

A convenience wrapper over odeint for the common case of a single interval with steps uniform steps; differentiable in y0 and theta.

estimate_dynamics ¶

estimate_dynamics(
    dynamics: Callable[[Array, Any, Any], Any],
    y0: Any,
    ts: Array,
    data: Array,
    theta0: Any,
    *,
    observe: Callable[[Any], Array] | None = None,
    weights: Array | None = None,
    integrator: str = "rk4",
    substeps: int = 2,
    max_iter: int = 100,
) -> DynamicEstimateResult

Fit parameters theta so the simulated trajectory matches data.

Minimizes sum (w (observe(traj) - data))^2 by Levenberg-Marquardt, with the trajectory produced by integrating dynamics(t, y, theta) from y0 over ts. observe maps the trajectory pytree to the measured quantity (defaults to the trajectory itself); weights optionally scales residuals.

Returns:

Type	Description
`DynamicEstimateResult`	A `DynamicEstimateResult`.

optimal_control ¶

optimal_control(
    dynamics: Callable[[Array, Any, Any, Any], Any],
    y0: Any,
    ts: Array,
    u_init: Array,
    stage_cost: Callable[[Array, Any, Any, Any], Array],
    *,
    terminal_cost: Callable[[Any, Any], Array]
    | None = None,
    theta: Any = None,
    bounds: tuple[Any, Any] | None = None,
    method: str = "bfgs",
    integrator: str = "rk4",
    substeps: int = 2,
    max_iter: int = 100,
) -> OptimalControlResult

Find the input trajectory minimizing a running-plus-terminal cost.

The control is piecewise constant on the ts grid: u_init has one entry (scalar or vector) per interval. The augmented state (y, J) integrates the running cost J' = stage_cost(t, y, u, theta) alongside the dynamics y' = dynamics(t, y, u, theta), and the total J(t_f) + terminal_cost(y_f) is minimized over the control with gradients straight through the simulation.

Parameters:

Name	Type	Description	Default
`dynamics`	`Callable[[Array, Any, Any, Any], Any]`	`dynamics(t, y, u, theta) -> dy` (note the explicit input `u`).	required
`y0`	`Any`	Initial state pytree.	required
`ts`	`Array`	Output/decision grid (length `N`); the control has `N - 1` intervals.	required
`u_init`	`Array`	Initial control guess, shape `(N - 1,)` or `(N - 1, k)`.	required
`stage_cost`	`Callable[[Array, Any, Any, Any], Array]`	Running cost `stage_cost(t, y, u, theta) -> ()`.	required
`terminal_cost`	`Callable[[Any, Any], Array] \| None`	Optional terminal cost `terminal_cost(y_f, theta) -> ()`.	`None`
`theta`	`Any`	Optional fixed parameter pytree.	`None`
`bounds`	`tuple[Any, Any] \| None`	Optional `(lower, upper)` box on the control.	`None`
`method`	`str`	Unconstrained optimizer used over the control (e.g. `"bfgs"`).	`'bfgs'`
`integrator`	`str`	ODE integration scheme for the augmented state (e.g. `"rk4"`).	`'rk4'`
`substeps`	`int`	Integration substeps taken between consecutive `ts` points.	`2`
`max_iter`	`int`	Maximum number of optimizer iterations.	`100`

Returns:

Type	Description
`OptimalControlResult`	An `OptimalControlResult`.

tune_pid ¶

tune_pid(
    response: Callable[[Any], Array],
    gains0: Any,
    setpoint: ArrayLike,
    ts: Array,
    *,
    objective: str = "iae",
    bounds: tuple[Any, Any] | None = None,
    method: str = "bfgs",
    max_iter: int = 100,
) -> OptimizeResult

Tune controller gains by minimizing a closed-loop error integral.

response(gains) must build the closed loop from the gains pytree, simulate it, and return the controlled-variable trajectory sampled on ts. This function then minimizes the chosen error integral ("iae", "ise" or "itae") against setpoint over the gains: gradients flow through the whole simulated loop, so the tune is exact first-order, not a grid search.

Parameters:

Name	Type	Description	Default
`response`	`Callable[[Any], Array]`	`response(gains) -> pv_trajectory` (length `len(ts)`).	required
`gains0`	`Any`	Initial gains pytree (e.g. a dict `{"kc": ..., "tau_i": ...}`).	required
`setpoint`	`ArrayLike`	Target value for the controlled variable.	required
`ts`	`Array`	Time grid matching `response`'s output.	required
`objective`	`str`	Error integral to minimize.	`'iae'`
`bounds`	`tuple[Any, Any] \| None`	Optional `(lower, upper)` box on the gains (recommended: keeps gains positive).	`None`
`method`	`str`	Unconstrained optimizer used over the gains (e.g. `"bfgs"`).	`'bfgs'`
`max_iter`	`int`	Maximum number of optimizer iterations.	`100`

Returns:

Type	Description
`OptimizeResult`	The `fugacio.sim.OptimizeResult`; `result.x` is the tuned gains.

Control¶

control ¶

Classical and differentiable process control for Fugacio.

This subpackage supplies the control side of dynamic simulation:

pid: a realizable, anti-windup PID controller whose gains are a differentiable pytree (so they can be tuned by gradient descent), carried as ODE states inside a dynamic flowsheet;
blocks: linear blocks (first/second order, FOPDT, lead-lag) with both analytic step responses and state-space realizations, plus static actuator nonlinearities;
metrics: step-response performance metrics (overshoot, rise/settling time, IAE/ISE/ITAE) for reporting and as tuning objectives;
tuning: FOPDT model identification and the classical tuning rules (Ziegler-Nichols, Cohen-Coon, IMC/lambda, AMIGO);
linearize: exact autodiff linearization of a nonlinear plant into a StateSpace, with poles, Bode response, and controllability/observability.

Modules:

Name	Description
`blocks`	Linear dynamic blocks and analytic step responses for control studies.
`linearize`	Linearize a nonlinear dynamic model about an operating point, by autodiff.
`metrics`	Time-domain performance metrics for a response trajectory.
`pid`	A differentiable PID controller in realizable, anti-windup form.
`tuning`	PID tuning rules and first-order-plus-dead-time model identification.

Classes:

Name	Description
`StateSpace`	A linear time-invariant state-space model `x' = A x + B u`, `y = C x + D u`.
`StepInfo`	A bundle of step-response metrics.
`PID`	A filtered, anti-windup PID controller (a differentiable pytree).
`PIDState`	Internal controller state carried through the dynamic simulation.
`FOPDTModel`	A first-order-plus-dead-time process model `K e^{-L s} / (tau s + 1)`.

Functions:

Name	Description
`dead_band`	Symmetric dead band of total width `2*width` centred at zero.
`first_order_ss`	State-space `(A, B, C, D)` of `K/(tau s + 1)` (one state).
`first_order_step`	Response of `K/(tau s + 1)` to a step of size `u` at `t = 0`.
`fopdt_step`	Response of a first-order-plus-dead-time process to a step of size `u`.
`lead_lag`	Step response of a lead-lag compensator `K (tau_lead s + 1)/(tau_lag s + 1)`.
`rate_limit`	Limit a desired rate of change to `+/- max_rate`.
`saturate`	Clamp `u` to `[u_min, u_max]`.
`second_order_ss`	State-space `(A, B, C, D)` of the standard second-order block (two states).
`second_order_step`	Step response of `K wn^2 / (s^2 + 2 zeta wn s + wn^2)` (size `u`).
`bode`	Bode magnitude (dB) and phase (degrees) of a SISO system over `omega`.
`controllability`	Controllability matrix `[B, A B, ..., A^{n-1} B]`.
`dc_gain`	Steady-state gain `-C A^{-1} B + D` (the response to a unit step at `t -> inf`).
`frequency_response`	Complex frequency response `H(j omega) = C (j omega I - A)^{-1} B + D`.
`is_controllable`	Whether the controllability matrix has full state rank.
`is_observable`	Whether the observability matrix has full state rank.
`is_stable`	Whether every pole has strictly negative real part (asymptotic stability).
`observability`	Observability matrix `[C; C A; ...; C A^{n-1}]`.
`poles`	Eigenvalues of `A` (the system poles).
`iae`	Integral of the absolute error, `integral \|sp - y\| dt`.
`ise`	Integral of the squared error, `integral (sp - y)^2 dt`.
`itae`	Integral of time-weighted absolute error, `integral t \|sp - y\| dt`.
`overshoot`	Fractional overshoot of the peak beyond the setpoint (0 if none).
`peak_time`	Time of the response extremum in the step direction.
`rise_time`	Time to rise from `lo` to `hi` fraction of the step (`nan` if never reached).
`settling_time`	Last time the response leaves the `+/- tol` band around the setpoint.
`steady_state_error`	Offset of the final value from the setpoint, `sp - y[-1]`.
`step_info`	Compute the standard step-response metrics for a trajectory `(t, y)`.
`p_only`	Convenience constructor for a proportional-only controller.
`pi`	Convenience constructor for a PI controller (no derivative action).
`amigo`	AMIGO (Astrom-Hagglund) robust tuning from an FOPDT model.
`cohen_coon`	Cohen-Coon tuning from an FOPDT model (better than ZN for larger `L/tau`).
`fit_fopdt`	Identify an FOPDT model from a step response `(t, y)` to an input step `u_step`.
`imc_tuning`	Internal-model-control (lambda) tuning from an FOPDT model.
`ziegler_nichols`	Ziegler-Nichols open-loop (reaction-curve) tuning from an FOPDT model.

StateSpace ¶

Bases: NamedTuple

A linear time-invariant state-space model x' = A x + B u, y = C x + D u.

Attributes:

Name	Type	Description
`a`	`Array`	State matrix `(n, n)`.
`b`	`Array`	Input matrix `(n, m)`.
`c`	`Array`	Output matrix `(p, n)`.
`d`	`Array`	Feedthrough matrix `(p, m)`.

StepInfo ¶

Bases: NamedTuple

A bundle of step-response metrics.

Attributes:

Name	Type	Description
`overshoot`	`Array`	Fractional overshoot beyond the setpoint.
`peak_time`	`Array`	Time of the response extremum.
`rise_time`	`Array`	10-90% rise time.
`settling_time`	`Array`	2% settling time.
`steady_state_error`	`Array`	Final offset `sp - y[-1]`.
`iae`	`Array`	Integral of absolute error.

PID `dataclass` ¶

PID(
    kc: ArrayLike,
    tau_i: ArrayLike = inf,
    tau_d: ArrayLike = 0.0,
    beta: ArrayLike = 1.0,
    gamma: ArrayLike = 0.0,
    u_bias: ArrayLike = 0.0,
    u_min: ArrayLike = -inf,
    u_max: ArrayLike = inf,
    tau_t: ArrayLike = 0.0,
    n_filter: float = 10.0,
    direction: str = "reverse",
)

A filtered, anti-windup PID controller (a differentiable pytree).

Attributes:

Name	Type	Description
`kc`	`ArrayLike`	Proportional gain (output units per measurement unit).
`tau_i`	`ArrayLike`	Integral (reset) time, s. Use `inf` to disable integral action.
`tau_d`	`ArrayLike`	Derivative (rate) time, s. `0` disables derivative action.
`beta`	`ArrayLike`	Setpoint weight on the proportional term (`0..1`; `1` is classic).
`gamma`	`ArrayLike`	Setpoint weight on the derivative term (usually `0`, derivative on measurement only, to avoid derivative kick).
`u_bias`	`ArrayLike`	Output bias / nominal output at zero error (manual reset).
`u_min`	`ArrayLike`	Lower output limit (saturation).
`u_max`	`ArrayLike`	Upper output limit (saturation).
`tau_t`	`ArrayLike`	Anti-windup tracking time, s. `<= 0` selects an automatic value (`sqrt(tau_i * tau_d)` when derivative is active, else `tau_i`).
`n_filter`	`float`	Derivative-filter ratio (filter time is `tau_d / n_filter`).
`direction`	`str`	`"reverse"` (default) or `"direct"` acting; `"direct"` flips the error sign.

Methods:

Name	Description
`init_state`	Initial controller state for bumpless start at output `u0` (default bias).
`output`	Saturated controller output for the current state, setpoint and measurement.
`derivative`	Time derivative of the controller state (for embedding in a flowsheet ODE).
`step`	Discrete one-step update by explicit Euler; returns `(output, new_state)`.

init_state ¶

init_state(
    pv0: ArrayLike, u0: ArrayLike | None = None
) -> PIDState

Initial controller state for bumpless start at output u0 (default bias).

The integral term is preloaded so the controller's initial output equals u0 (or u_bias if u0 is None) when the measurement is at pv0 and the setpoint equals pv0; the derivative filter starts at the measurement so the initial derivative action is zero.

output ¶

output(
    state: PIDState, setpoint: ArrayLike, pv: ArrayLike
) -> Array

Saturated controller output for the current state, setpoint and measurement.

derivative ¶

derivative(
    state: PIDState, setpoint: ArrayLike, pv: ArrayLike
) -> PIDState

Time derivative of the controller state (for embedding in a flowsheet ODE).

Returns (di/dt, dx_d/dt) as a PIDState. The integral derivative includes the back-calculation anti-windup term so it stops winding up while the output is saturated.

step ¶

step(
    state: PIDState,
    setpoint: ArrayLike,
    pv: ArrayLike,
    dt: ArrayLike,
) -> tuple[Array, PIDState]

Discrete one-step update by explicit Euler; returns (output, new_state).

For standalone digital-controller loops. Inside a continuous dynamic simulation prefer derivative so the controller integrates with the same solver as the plant.

PIDState ¶

Bases: NamedTuple

Internal controller state carried through the dynamic simulation.

Attributes:

Name	Type	Description
`i`	`Array`	Integral term contribution to the output (already scaled by the integral gain), in output units.
`x_d`	`Array`	First-order-filtered measurement used to form the derivative term.

FOPDTModel ¶

Bases: NamedTuple

A first-order-plus-dead-time process model K e^{-L s} / (tau s + 1).

Attributes:

Name	Type	Description
`gain`	`Array`	Steady-state gain `K` (output units per input unit).
`tau`	`Array`	Time constant `tau` (s).
`dead_time`	`Array`	Dead time / transport delay `L` (s).

dead_band ¶

dead_band(e: ArrayLike, width: ArrayLike) -> Array

Symmetric dead band of total width 2*width centred at zero.

first_order_ss ¶

first_order_ss(
    gain: ArrayLike, tau: ArrayLike
) -> tuple[Array, Array, Array, Array]

State-space (A, B, C, D) of K/(tau s + 1) (one state).

first_order_step ¶

first_order_step(
    t: ArrayLike,
    gain: ArrayLike,
    tau: ArrayLike,
    *,
    u: ArrayLike = 1.0,
) -> Array

Response of K/(tau s + 1) to a step of size u at t = 0.

y(t) = K u (1 - exp(-t/tau)) for t >= 0.

fopdt_step ¶

fopdt_step(
    t: ArrayLike,
    gain: ArrayLike,
    tau: ArrayLike,
    dead_time: ArrayLike,
    *,
    u: ArrayLike = 1.0,
) -> Array

Response of a first-order-plus-dead-time process to a step of size u.

y(t) = K u (1 - exp(-(t - L)/tau)) for t >= L (the dead time L), else 0. The workhorse model for PID tuning.

lead_lag ¶

lead_lag(
    t: ArrayLike,
    gain: ArrayLike,
    tau_lead: ArrayLike,
    tau_lag: ArrayLike,
    *,
    u: ArrayLike = 1.0,
) -> Array

Step response of a lead-lag compensator K (tau_lead s + 1)/(tau_lag s + 1).

rate_limit ¶

rate_limit(
    du_desired: ArrayLike, max_rate: ArrayLike
) -> Array

Limit a desired rate of change to +/- max_rate.

saturate ¶

saturate(
    u: ArrayLike, u_min: ArrayLike, u_max: ArrayLike
) -> Array

Clamp u to [u_min, u_max].

second_order_ss ¶

second_order_ss(
    gain: ArrayLike, wn: ArrayLike, zeta: ArrayLike
) -> tuple[Array, Array, Array, Array]

State-space (A, B, C, D) of the standard second-order block (two states).

second_order_step ¶

second_order_step(
    t: ArrayLike,
    gain: ArrayLike,
    wn: ArrayLike,
    zeta: ArrayLike,
    *,
    u: ArrayLike = 1.0,
) -> Array

Step response of K wn^2 / (s^2 + 2 zeta wn s + wn^2) (size u).

Handles the under-, critically-, and over-damped regimes with a single smooth expression (the under/over branches are selected by jax.numpy.where, so the result is differentiable in zeta through zeta = 1 as well).

bode ¶

bode(ss: StateSpace, omega: Array) -> tuple[Array, Array]

Bode magnitude (dB) and phase (degrees) of a SISO system over omega.

controllability ¶

controllability(ss: StateSpace) -> Array

Controllability matrix [B, A B, ..., A^{n-1} B].

dc_gain ¶

dc_gain(ss: StateSpace) -> Array

Steady-state gain -C A^{-1} B + D (the response to a unit step at t -> inf).

frequency_response ¶

frequency_response(ss: StateSpace, omega: Array) -> Array

Complex frequency response H(j omega) = C (j omega I - A)^{-1} B + D.

Returns an array of shape (len(omega), p, m) of complex transfer-function values.

is_controllable ¶

is_controllable(ss: StateSpace) -> Array

Whether the controllability matrix has full state rank.

is_observable ¶

is_observable(ss: StateSpace) -> Array

Whether the observability matrix has full state rank.

is_stable ¶

is_stable(ss: StateSpace) -> Array

Whether every pole has strictly negative real part (asymptotic stability).

observability ¶

observability(ss: StateSpace) -> Array

Observability matrix [C; C A; ...; C A^{n-1}].

poles ¶

poles(ss: StateSpace) -> Array

Eigenvalues of A (the system poles).

iae ¶

iae(t: Array, y: Array, setpoint: ArrayLike) -> Array

Integral of the absolute error, integral |sp - y| dt.

ise ¶

ise(t: Array, y: Array, setpoint: ArrayLike) -> Array

Integral of the squared error, integral (sp - y)^2 dt.

itae ¶

itae(t: Array, y: Array, setpoint: ArrayLike) -> Array

Integral of time-weighted absolute error, integral t |sp - y| dt.

overshoot ¶

overshoot(
    y: Array,
    setpoint: ArrayLike,
    *,
    y0: ArrayLike | None = None,
) -> Array

Fractional overshoot of the peak beyond the setpoint (0 if none).

(y_peak - sp) / (sp - y0) for a positive step; y0 defaults to the first sample. Returns 0 when the response does not exceed the setpoint.

peak_time ¶

peak_time(t: Array, y: Array, setpoint: ArrayLike) -> Array

Time of the response extremum in the step direction.

rise_time ¶

rise_time(
    t: Array,
    y: Array,
    setpoint: ArrayLike,
    *,
    lo: float = 0.1,
    hi: float = 0.9,
) -> Array

Time to rise from lo to hi fraction of the step (nan if never reached).

settling_time ¶

settling_time(
    t: Array,
    y: Array,
    setpoint: ArrayLike,
    *,
    tol: float = 0.02,
) -> Array

Last time the response leaves the +/- tol band around the setpoint.

Returns the time after which |y - sp| <= tol * |sp - y0| holds for the rest of the record (t[-1] if it never settles within the record).

steady_state_error ¶

steady_state_error(y: Array, setpoint: ArrayLike) -> Array

Offset of the final value from the setpoint, sp - y[-1].

step_info ¶

step_info(
    t: Array,
    y: Array,
    setpoint: ArrayLike,
    *,
    settle_tol: float = 0.02,
) -> StepInfo

Compute the standard step-response metrics for a trajectory (t, y).

p_only ¶

p_only(kc: ArrayLike, **kwargs: ArrayLike) -> PID

Convenience constructor for a proportional-only controller.

pi ¶

pi(
    kc: ArrayLike, tau_i: ArrayLike, **kwargs: ArrayLike
) -> PID

Convenience constructor for a PI controller (no derivative action).

amigo ¶

amigo(
    model: FOPDTModel,
    *,
    controller: str = "PI",
    **kwargs: ArrayLike,
) -> PID

AMIGO (Astrom-Hagglund) robust tuning from an FOPDT model.

cohen_coon ¶

cohen_coon(
    model: FOPDTModel,
    *,
    controller: str = "PID",
    **kwargs: ArrayLike,
) -> PID

Cohen-Coon tuning from an FOPDT model (better than ZN for larger L/tau).

fit_fopdt ¶

fit_fopdt(
    t: Array,
    y: Array,
    u_step: ArrayLike = 1.0,
    *,
    guess: FOPDTModel | None = None,
) -> FOPDTModel

Identify an FOPDT model from a step response (t, y) to an input step u_step.

Fits (K, tau, L) by Levenberg-Marquardt least squares against the analytic FOPDT step response. A data-driven initial guess is used when guess is omitted (gain from the final value; tau and L from the 28.3% / 63.2% response times). The fit is differentiable in the data.

imc_tuning ¶

imc_tuning(
    model: FOPDTModel,
    *,
    tau_c: ArrayLike | None = None,
    controller: str = "PI",
    **kwargs: ArrayLike,
) -> PID

Internal-model-control (lambda) tuning from an FOPDT model.

tau_c is the desired closed-loop time constant; it defaults to max(0.1 tau, 0.8 L), a robust middle-of-the-road choice. Larger tau_c gives a slower but more robust loop.

ziegler_nichols ¶

ziegler_nichols(
    model: FOPDTModel,
    *,
    controller: str = "PID",
    **kwargs: ArrayLike,
) -> PID

Ziegler-Nichols open-loop (reaction-curve) tuning from an FOPDT model.

controller is "P", "PI" or "PID". Aggressive (quarter-amplitude decay) tuning, a classic baseline, often detuned in practice.

Dynamics & process control¶

Dynamics¶

dynamics ¶

DynamicFlowsheet dataclass ¶

feed ¶

add ¶

control ¶

manipulate ¶

initial_state ¶

rhs ¶

simulate ¶

DynamicResult ¶

measurement ¶

ODEResult ¶

DynamicEstimateResult ¶

OptimalControlResult ¶

DynamicCSTR dataclass ¶

DynamicFlash dataclass ¶

DynamicUnit ¶

initial_state abstractmethod ¶

evaluate abstractmethod ¶

measured abstractmethod ¶

GasReceiver dataclass ¶

LevelTank dataclass ¶

MixingTank dataclass ¶

ThermalMass dataclass ¶

simulate ¶

odeint ¶

odeint_final ¶

estimate_dynamics ¶

optimal_control ¶

tune_pid ¶

Control¶

control ¶

StateSpace ¶

StepInfo ¶

PID dataclass ¶

init_state ¶

output ¶

derivative ¶

step ¶

PIDState ¶

FOPDTModel ¶

dead_band ¶

first_order_ss ¶

first_order_step ¶

fopdt_step ¶

lead_lag ¶

rate_limit ¶

saturate ¶

second_order_ss ¶

second_order_step ¶

bode ¶

controllability ¶

dc_gain ¶

frequency_response ¶

is_controllable ¶

is_observable ¶

is_stable ¶

observability ¶

poles ¶

iae ¶

ise ¶

itae ¶

overshoot ¶

peak_time ¶

rise_time ¶

settling_time ¶

steady_state_error ¶

step_info ¶

p_only ¶

pi ¶

amigo ¶

cohen_coon ¶

fit_fopdt ¶

imc_tuning ¶

ziegler_nichols ¶

DynamicFlowsheet `dataclass` ¶

DynamicCSTR `dataclass` ¶

DynamicFlash `dataclass` ¶

initial_state `abstractmethod` ¶

evaluate `abstractmethod` ¶

measured `abstractmethod` ¶

GasReceiver `dataclass` ¶

LevelTank `dataclass` ¶

MixingTank `dataclass` ¶

ThermalMass `dataclass` ¶

PID `dataclass` ¶