Anatomy of a Symbolic System

Stochastic Batch Reactors

Author

Gil Benezer

Published

January 16, 2026

This tutorial introduces the anatomy of a symbolic control system using a stochastic batch reactor as a worked example. A batch reactor is a closed chemical system where reactants $A$ and $B$ undergo sequential reactions governed by Arrhenius kinetics, with temperature controlled via external heating $Q$.

System Overview

The system demonstrates key cdesym concepts:

State variables: Concentrations $(C_A, C_B)$ and temperature $(T)$
Control input: Heat input $(Q)$
Drift dynamics: Deterministic reaction rates and heat transfer
Diffusion dynamics: Additive noise on each state, representing measurement or process uncertainty

The reaction sequence $A \to B \to C$ follows first-order Arrhenius kinetics:

\[ \begin{aligned} r_1 &= k_1 \cdot C_A \cdot \exp(\frac{-E_1}{T}) \\ r_2 &= k_2 \cdot C_B \cdot \exp(\frac{-E_2}{T}) \end{aligned} \]

The deterministic dynamics are:

\[ \begin{aligned} \frac{dC_A}{dt} &= -r_1 \\ \frac{dC_B}{dt} &= r_1 - r_2 \\ \frac{dT}{dt} &= Q - \alpha(T - T_{\text{amb}}) \end{aligned} \]

Itô vs Stratonovich SDEs

The stochastic formulation in this example uses an Itô SDE of the form:

\[dX_t = f(X_t, u_t)\,dt + g(X_t)\,dW_t\]

where $W_t$ is a Wiener process (Brownian motion). The key distinction from Stratonovich SDEs lies in how the stochastic integral is defined:

Property	Itô	Stratonovich
Evaluation point	Left endpoint of interval	Midpoint of interval
Martingale property	Yes (integral is a martingale)	No
Chain rule	Requires Itô’s lemma with correction term	Standard chain rule applies
Typical use	Finance, control theory, probability	Physics, systems with colored noise limits

Itô’s lemma for a function $Y_t = h(X_t)$ includes an extra term:

\[dY_t = h'(X_t)\,dX_t + \frac{1}{2}h''(X_t)\,g(X_t)^2\,dt\]

The second term (absent in ordinary calculus) arises because Brownian paths have infinite quadratic variation.

When the distinction matters: For state-dependent (multiplicative) noise like $g(X_t) = \sigma X_t$, the two interpretations give different results. For additive noise where $g$ is constant, they coincide. This tutorial uses additive noise, so the choice is moot here—but understanding the difference is essential when modeling state-dependent noise.

Noise Structure

The stochastic batch reactor extends the deterministic model with additive process noise:

\[ \begin{aligned} dC_A &= (-r_1)\,dt + \sigma_A\,dW_A \\ dC_B &= (r_1 - r_2)\,dt + \sigma_B\,dW_B \\ dT &= (Q - \alpha(T - T_{\text{amb}}))\,dt + \sigma_T\,dW_T \end{aligned} \]

Three independent Wiener processes model distinct physical noise sources: feed variability ($\sigma_A$, $\sigma_B$) and heat transfer fluctuations ($\sigma_T$). Trajectories fluctuate around the deterministic equilibrium with standard deviation growing as $\sqrt{t}$.

What This Tutorial Covers

We define two equilibria:

Complete conversion (stable): All reactants consumed, temperature at ambient
Initial setpoint (unstable): Starting conditions for batch operation

The tutorial walks through defining _f_sym (drift), diffusion_expr (noise structure), how to define optional _h_sym (output map) and output variables if needed, then covers equilibrium setup and trajectory visualization.

System Definition

The full system definition is as follows:

Code

from typing import Optional

import numpy as np
import sympy as sp

from cdesym import ContinuousStochasticSystem

class ContinuousStochasticBatchReactor(ContinuousStochasticSystem):
    def define_system(
        self,
        k1_val: float = 0.5,
        k2_val: float = 0.3,
        E1_val: float = 1000.0,
        E2_val: float = 1500.0,
        alpha_val: float = 0.1,
        T_amb_val: float = 300.0,
        sigma_A: float = 0.01,
        sigma_B: float = 0.01,
        sigma_T: float = 1.0,
        C_A0: Optional[float] = None,
        T0: Optional[float] = None,
    ):
        # Store initial conditions and parameter values for later use
        # (See "What to Store as Attributes" callout below for detailed guidelines)
        self.C_A0 = C_A0
        self.T0 = T0
        self.alpha_val = alpha_val
        self.T_amb_val = T_amb_val

        # State and control variables (local - not stored as attributes)
        C_A, C_B, T = sp.symbols("C_A C_B T", real=True, positive=True)
        Q = sp.symbols("Q", real=True)

        # Deterministic Parameters (kinetics and heat transfer)
        # These are local variables - used as keys in self.parameters dict
        k1, k2, E1, E2, alpha, T_amb = sp.symbols(
            "k1 k2 E1 E2 alpha T_amb",
            real=True,
            positive=True,
        )

        # Stochastic Parameters (noise intensities)
        sigma_A_sym = sp.symbols("sigma_A", real=True, positive=True)
        sigma_B_sym = sp.symbols("sigma_B", real=True, positive=True)
        sigma_T_sym = sp.symbols("sigma_T", real=True, positive=True)

        # Parameters dictionary: Symbol → Value mapping
        # Keys are the symbol objects created above
        self.parameters = {
            k1: k1_val,
            k2: k2_val,
            E1: E1_val,
            E2: E2_val,
            alpha: alpha_val,
            T_amb: T_amb_val,
            sigma_A_sym: sigma_A,
            sigma_B_sym: sigma_B,
            sigma_T_sym: sigma_T,
        }

        self.state_vars = [C_A, C_B, T]
        self.control_vars = [Q]
        self.output_vars = []
        self.order = 1

        # Reaction rates (Arrhenius kinetics)
        r1 = k1 * C_A * sp.exp(-E1 / T)
        r2 = k2 * C_B * sp.exp(-E2 / T)

        # DRIFT (Deterministic part - same as deterministic reactor)
        dC_A_dt = -r1
        dC_B_dt = r1 - r2
        dT_dt = Q - alpha * (T - T_amb)

        self._f_sym = sp.Matrix([dC_A_dt, dC_B_dt, dT_dt])

        # DIFFUSION (Stochastic part - additive noise)
        # Diagonal matrix: three independent Wiener processes
        self.diffusion_expr = sp.Matrix(
            [
                [sigma_A_sym, 0, 0],
                [0, sigma_B_sym, 0],
                [0, 0, sigma_T_sym],
            ],
        )

        # Itô SDE interpretation
        self.sde_type = "ito"

        # Output: Full state measurement (with potential noise in practice)
        # technically optional
        self._h_sym = sp.Matrix([C_A, C_B, T])

    def setup_equilibria(self):
        # Use stored parameter values directly
        # (Don't try to look up via sp.symbols() - that creates NEW symbol objects!)
        T_amb = self.T_amb_val

        # Complete conversion equilibrium (deterministic part)
        self.add_equilibrium(
            "complete",
            x_eq=np.array([0.0, 0.0, T_amb]),
            u_eq=np.array([0.0]),
            verify=True,
            stability="stable",
            notes="Equilibrium of deterministic part (drift). Stochastic trajectories "
            "fluctuate around this point with variance growing over time.",
        )

        # Initial condition (if provided)
        if self.C_A0 is not None and self.T0 is not None:
            alpha = self.alpha_val
            Q_init = alpha * (self.T0 - T_amb)

            self.add_equilibrium(
                "initial",
                x_eq=np.array([self.C_A0, 0.0, self.T0]),
                u_eq=np.array([Q_init]),
                verify=False,
                stability="unstable",
                notes="Initial state setpoint (deterministic part). Stochastic trajectories "
                "will fluctuate with std ~ [σ_A·√t, σ_B·√t, σ_T·√t].",
            )
            self.set_default_equilibrium("initial")
        else:
            self.set_default_equilibrium("complete")

The first step in system construction is to subclass the proper symbolic system and create the define_system method signature.

class ContinuousStochasticBatchReactor(ContinuousStochasticSystem):
    def define_system(
        self,
        k1_val: float = 0.5,
        k2_val: float = 0.3,
        E1_val: float = 1000.0,
        E2_val: float = 1500.0,
        alpha_val: float = 0.1,
        T_amb_val: float = 300.0,
        sigma_A: float = 0.01,
        sigma_B: float = 0.01,
        sigma_T: float = 1.0,
        C_A0: Optional[float] = None,
        T0: Optional[float] = None,
    ):

More often than not, the only arguments to this function will be the self argument and parameter values. It is not necessary to add type hints or default values, but it is recommended for debug purposes.

The next step in system definition is to define the (deterministic) symbolic variables to be used and place them in the proper fields. The main limitation currently is that it is not possible to explicitly define a symbolic variable for time. If it is heavily requested, we can implement it as a future feature.

# State and control variables
C_A, C_B, T = sp.symbols("C_A C_B T", real=True, positive=True)
Q = sp.symbols("Q", real=True)

# Deterministic Parameters (kinetics and heat transfer)
k1, k2, E1, E2, alpha, T_amb = sp.symbols(
    "k1 k2 E1 E2 alpha T_amb",
    real=True,
    positive=True,
)

# Stochastic Parameters (noise intensities)
sigma_A_sym = sp.symbols("sigma_A", real=True, positive=True)
sigma_B_sym = sp.symbols("sigma_B", real=True, positive=True)
sigma_T_sym = sp.symbols("sigma_T", real=True, positive=True)

self.parameters = {
    k1: k1_val,
    k2: k2_val,
    E1: E1_val,
    E2: E2_val,
    alpha: alpha_val,
    T_amb: T_amb_val,
    sigma_A_sym: sigma_A,
    sigma_B_sym: sigma_B,
    sigma_T_sym: sigma_T,
}

self.state_vars = [C_A, C_B, T]
self.control_vars = [Q]
self.output_vars = []

The necessary symbolic variables include:

State variables
Control variables
Parameter variables

What to Store as Attributes vs Local Variables

A common source of confusion: when should you store something as self.something vs keep it as a local variable?

SymPy Symbols: Usually Local, Sometimes Stored

State and control symbols → Keep LOCAL (don’t store as attributes):

# CORRECT: Local variables only
C_A, C_B, T = sp.symbols("C_A C_B T", real=True)
Q = sp.symbols("Q", real=True)

# WRONG: Don't store these as attributes
# self.C_A_sym = sp.symbols("C_A")  # Unnecessary!

Why keep local? These symbols are only used to build self._f_sym and self._h_sym. Once those are built, the framework tracks symbols via self.state_vars and self.control_vars lists. Storing individual symbols is redundant.

Parameter symbols → Keep LOCAL (unless using Solution 2 for equilibria):

# Solution 1 (tutorial approach): Local symbols
T_amb = sp.symbols('T_amb', positive=True)
self.parameters = {T_amb: T_amb_val}  # Symbol used as key only

# Solution 2 (production code): Store symbol references
self.T_amb_sym = sp.symbols('T_amb', positive=True)
self.parameters = {self.T_amb_sym: T_amb_val}  # Can look up later

Parameter Values: Store When Needed Later

Store parameter values if used in setup_equilibria() or custom methods:

def define_system(self, T_amb_val=300.0, alpha_val=0.1, ...):
    # Store values needed in setup_equilibria
    self.T_amb_val = T_amb_val
    self.alpha_val = alpha_val
    
    # Create symbols and populate parameters dict
    T_amb = sp.symbols('T_amb', positive=True)
    self.parameters = {T_amb: T_amb_val, ...}

Don’t store if only used in symbolic expressions:

def define_system(self, k1_val=0.5, k2_val=0.3, ...):
    # These are only used in _f_sym, not in setup_equilibria
    # No need to store as attributes
    k1 = sp.symbols('k1', positive=True)
    k2 = sp.symbols('k2', positive=True)
    
    self.parameters = {k1: k1_val, k2: k2_val}
    self._f_sym = sp.Matrix([k1 * C_A * sp.exp(-E1/T)])  # Used here only

Custom Fields: Always Store

Store data needed for conditional logic or custom methods:

def define_system(self, C_A0=None, T0=None, ...):
    # Store for use in setup_equilibria
    self.C_A0 = C_A0
    self.T0 = T0
    
    # Store custom application data
    self.reaction_pathway = "A->B->C"
    self.reactor_volume = 10.0  # liters

Decision Tree

Should I store foo as self.foo?

Is it a SymPy symbol?
- State/control symbol? → No (use self.state_vars / self.control_vars lists)
- Parameter symbol? → No (unless using Solution 2 for parameter lookup)
Is it a parameter value?
- Used in setup_equilibria()? → Yes (self.param_val)
- Only used in _f_sym? → No (just put in self.parameters dict)
Is it a custom field?
- Used in logic or custom methods? → Yes (self.C_A0, self.reactor_volume)
- Temporary calculation? → No (keep as local variable)

Rule of thumb: If you’ll need it in setup_equilibria() or a custom method, store it as self.something. Otherwise, keep it local.

Output Variables and the Output Map

There are two related but distinct concepts for system outputs:

self.output_vars: A list of symbolic variables representing computed outputs—quantities derived from state variables but not states themselves. For example, if your states are position and velocity, an output variable might be kinetic energy. Do not add state variables to this list.
self._h_sym: The output map $y = h(x)$, a symbolic column vector defining what the system actually outputs. This can include:
- State variables (for full or partial state observation)
- Output variables (computed quantities)
- Any symbolic expression of states

Defaults: self.output_vars = [] (no computed outputs) and self._h_sym = state vector (full state observability).

Example: For a system with states $[x, v]$ and a computed output $E = \frac{1}{2}mv^2$:

self.output_vars = [E]
self._h_sym = sp.Matrix([x, E]) outputs position and energy (but not velocity)

Autonomous Systems

An autonomous system has no external control inputs—its dynamics depend only on the current state (and possibly time). Examples include free oscillators, population dynamics, and unforced chemical reactions.

Required Specification

To define an autonomous system, you must set self.control_vars = [] (an empty list). Do not use None or omit the field—this will cause errors during system initialization.

# CORRECT: Autonomous system
self.control_vars = []

# WRONG: Will cause errors
self.control_vars = None
# self.control_vars = ...  (omitting the field)

Dynamics for Autonomous Systems

For autonomous systems, self._f_sym should depend only on state variables and parameters, not control variables:

# Example: Damped harmonic oscillator (autonomous)
x, v = sp.symbols('x v', real=True)
k, b, m = sp.symbols('k b m', positive=True)

self.state_vars = [x, v]
self.control_vars = []  # No control
self.parameters = {k: 10.0, b: 0.5, m: 1.0}

# Dynamics: dx/dt = v, dv/dt = -(k/m)*x - (b/m)*v
self._f_sym = sp.Matrix([
    v,
    -(k/m)*x - (b/m)*v
])

Calling Conventions

When evaluating or integrating autonomous systems:

The u parameter in system(x, u) should be None or omitted
Integration accepts u=None or no control specification
The framework automatically handles zero control internally

# Evaluating dynamics
dxdt = system(x)              # Recommended
dxdt = system(x, u=None)      # Also valid
dxdt = system(x, np.array([])) # Works but unnecessary

# Integration
result = system.integrate(x0=x_init, t_span=(0, 10))
# No need to specify u - it's ignored for autonomous systems

Converting Between Controlled and Autonomous

You can model the same physical system in different ways:

Controlled oscillator (with external force):

self.state_vars = [x, v]
self.control_vars = [u]  # External force
self._f_sym = sp.Matrix([v, -(k/m)*x - (b/m)*v + u/m])

Free oscillator (autonomous, special case u=0):

self.state_vars = [x, v]
self.control_vars = []  # No control
self._f_sym = sp.Matrix([v, -(k/m)*x - (b/m)*v])

Both are valid modeling choices. Use controlled form when you plan to apply inputs at some point; use autonomous form for analyzing free dynamics or when control is not part of the problem. Systems with control inputs expect at least a function that output constant zero input.

Time-Varying vs Autonomous

Important distinction:

Autonomous = no control inputs (self.control_vars = [])
Time-invariant = dynamics don’t explicitly depend on time
Time-varying = dynamics explicitly depend on time

A system can be autonomous but time-varying (e.g., $\dot{x} = -\sin(t) \cdot x$) or controlled but time-invariant (e.g., $\dot{x} = -x + u$). Most systems in control theory are both controlled and time-invariant.

Practical Examples

Autonomous systems:

Lotka-Volterra predator-prey: $\dot{x} = \alpha x - \beta xy$, $\dot{y} = \delta xy - \gamma y$
Van der Pol oscillator: $\ddot{x} + \mu(x^2-1)\dot{x} + x = 0$
Free batch reactor (no heating): $\dot{C}_A = -k_1 C_A e^{-E_1/T}$, $\dot{T} = -\alpha(T - T_{amb})$

Controlled systems:

Inverted pendulum with torque: $\ddot{\theta} = (g/L)\sin\theta + u/(mL^2)$
Batch reactor with heating (this tutorial): $\dot{T} = Q - \alpha(T - T_{amb})$
Linear regulator: $\dot{x} = Ax + Bu$

Parameter Dictionary Keys Must Be Symbolic Variables

A common mistake is to set the dictionary keys for self.parameters as strings, but for proper symbolic substitution and code generation, the keys must be SymPy Symbolic variables.

The final required step in system definition is to relate symbolic variables for the definition of the system (and possibly observation) dynamics.

# Reaction rates (Arrhenius kinetics)
r1 = k1 * C_A * sp.exp(-E1 / T)
r2 = k2 * C_B * sp.exp(-E2 / T)

# DRIFT (Deterministic part - same as deterministic reactor)
dC_A_dt = -r1
dC_B_dt = r1 - r2
dT_dt = Q - alpha * (T - T_amb)

self._f_sym = sp.Matrix([dC_A_dt, dC_B_dt, dT_dt])

# DIFFUSION (Stochastic part - additive noise)
# Diagonal matrix: three independent Wiener processes
self.diffusion_expr = sp.Matrix(
    [
        [sigma_A_sym, 0, 0],
        [0, sigma_B_sym, 0],
        [0, 0, sigma_T_sym],
    ],
)

Diffusion Matrix Structure

The diffusion matrix $G$ has dimensions $(n_{\text{states}} \times n_{\text{Wiener}})$, where each column corresponds to an independent Wiener process and each row to a state variable. The stochastic term in the SDE is $G\,dW_t$ where $dW_t$ is a vector of independent Wiener increments.

Why a matrix, not a vector? A vector would only allow one noise source affecting all states with fixed relative magnitudes. The matrix form enables:

Diagonal matrices (as above): Each state has its own independent noise source. The $i$-th Wiener process affects only the $i$-th state.
Off-diagonal elements: Represent correlated noise—when a single physical noise source affects multiple states. For example, if temperature fluctuations also induced concentration measurement errors:
```
self.diffusion_expr = sp.Matrix([
    [sigma_A_sym, 0, sigma_AT],  # C_A affected by W_A and W_T
    [0, sigma_B_sym, 0],
    [0, 0, sigma_T_sym],
])
```
Non-square matrices: You can have fewer Wiener processes than states (correlated noise) or more (redundant noise sources for modeling flexibility).

Dimensions expected: For $n$ state variables and $m$ independent Wiener processes, diffusion_expr should be an $n \times m$ SymPy Matrix.

# Output: Full state measurement (with potential noise in practice)
# technically optional
self._h_sym = sp.Matrix([C_A, C_B, T])

Pure Diffusion Systems

A pure diffusion system (also called a pure noise process) has no deterministic drift—dynamics are driven entirely by stochastic forcing. Examples include Brownian motion, geometric random walks, and models of measurement noise.

Mathematical Form

For pure diffusion systems, the SDE has zero drift:

\[dX_t = 0 \cdot dt + G(X_t)\,dW_t = G(X_t)\,dW_t\]

This means the state evolution is governed entirely by the diffusion term. The deterministic dynamics $f(x,u) = 0$ vanish.

Required Specification

To define a pure diffusion system, you must set self._f_sym to a zero vector with the same dimension as the number of state variables. Use SymPy’s zeros() function or explicitly construct the zero matrix:

# CORRECT: Pure diffusion system (3 states)
self._f_sym = sp.zeros(3, 1)  # Recommended approach

# ALSO CORRECT: Explicit zero matrix
self._f_sym = sp.Matrix([0, 0, 0])

# ALSO CORRECT: Using symbols and expressions that equal zero
self._f_sym = sp.Matrix([x - x, y - y, z - z])  # Works but unnecessary

# WRONG: Omitting _f_sym entirely (required field)
# self._f_sym = None  # Will cause errors

Example: 3D Brownian Motion

A particle undergoing pure diffusion in 3D space:

class BrownianMotion3D(ContinuousStochasticSystem):
    def define_system(self, sigma_x=1.0, sigma_y=1.0, sigma_z=1.0):
        # Position coordinates
        x, y, z = sp.symbols('x y z', real=True)
        
        # Diffusion coefficients
        sx, sy, sz = sp.symbols('sigma_x sigma_y sigma_z', positive=True)
        
        self.state_vars = [x, y, z]
        self.control_vars = []  # Autonomous (no control)
        self.parameters = {sx: sigma_x, sy: sigma_y, sz: sigma_z}
        
        # ZERO DRIFT: No deterministic motion
        self._f_sym = sp.zeros(3, 1)
        
        # DIFFUSION: Independent noise in each direction
        self.diffusion_expr = sp.Matrix([
            [sx, 0, 0],
            [0, sy, 0],
            [0, 0, sz]
        ])
        
        self.sde_type = "ito"

When to Use Pure Diffusion

Use pure diffusion when modeling:

Brownian motion: Pure random walks in continuous time
Measurement noise: Sensors with zero mean noise and no drift
Diffusion-only processes: Heat diffusion at steady-state, particle diffusion
Null models: Testing whether observed dynamics are purely random
Noise characterization: Estimating noise covariance without modeling dynamics

Don’t use pure diffusion when:

System has deterministic dynamics (even if noisy)
Modeling physical systems with forces/flows (use drift + diffusion)
Equilibrium points matter (pure diffusion has no equilibria except possibly at boundaries)

Physical Interpretation

For additive noise (constant $G$): \[dX_t = G\,dW_t\]

The solution is simply $X_t = X_0 + G\,W_t$, where $W_t$ is Brownian motion. The state performs a random walk around its initial condition with variance growing linearly: $\text{Var}(X_t) = GG^T \cdot t$.

For multiplicative noise (state-dependent $G(X)$): \[dX_t = G(X_t)\,dW_t\]

The dynamics are more complex, and the SDE interpretation (Itô vs Stratonovich) becomes critical.

Integration Considerations

Pure diffusion systems require:

Stochastic integrators: Methods like Euler-Maruyama, Milstein, or higher-order SDE schemes
Fine time steps: To resolve the noise structure accurately
Multiple realizations: Statistical properties emerge from ensemble averaging

# Example: Simulating 3D Brownian motion
brownian = BrownianMotion3D(sigma_x=1.0, sigma_y=1.0, sigma_z=1.0)

# Initial position at origin
x0 = np.array([0.0, 0.0, 0.0])

# Integrate using Euler-Maruyama
result = brownian.integrate(
    x0=x0,
    t_span=(0, 10),
    method="euler_maruyama",
    dt=0.01,  # Fine time step for noise resolution
    seed=42   # Reproducible noise
)

# For statistics, run multiple trajectories
trajectories = [
    brownian.integrate(x0=x0, t_span=(0, 10), dt=0.01, seed=i)
    for i in range(1000)
]

Drift + Diffusion vs Pure Diffusion

Most physical systems have both drift and diffusion:

Drift + Diffusion (this tutorial’s reactor):

# Deterministic part (drift)
self._f_sym = sp.Matrix([dC_A_dt, dC_B_dt, dT_dt])

# Stochastic part (diffusion)  
self.diffusion_expr = sp.Matrix([[σ_A, 0, 0], ...])

Pure Diffusion (rare in practice):

# No drift
self._f_sym = sp.zeros(3, 1)

# Only diffusion
self.diffusion_expr = sp.Matrix([[σ_A, 0, 0], ...])

The batch reactor in this tutorial has drift + diffusion because chemical reactions (drift) occur alongside process noise (diffusion).

Practical Applications

Pure Brownian Motion:

Modeling colloidal particle motion in fluids
Financial models with zero expected return
Null hypothesis testing in time series

State-Dependent Diffusion:

Volatility modeling in finance (stochastic volatility models)
Birth-death processes at steady-state
Fluctuation-dissipation in statistical mechanics

General Recommendation: If you’re unsure whether your system should be pure diffusion, it probably shouldn’t be. Most physical systems have deterministic dynamics (drift) plus noise (diffusion).

Required vs Optional System Definition Elements

Understanding what’s required versus optional helps you define systems efficiently while maintaining clarity. Here’s the complete breakdown:

Required Elements

Every define_system() method must define:

self.state_vars: List of state variable symbols
```
self.state_vars = [C_A, C_B, T]
```
self.control_vars: List of control variable symbols (empty list [] for autonomous systems—see the detailed Autonomous Systems callout)
```
self.control_vars = [Q]  # or [] for autonomous
```
self.parameters: Dictionary mapping symbolic parameters to numerical values
```
self.parameters = {k1: 0.5, k2: 0.3, ...}
```
self._f_sym: Symbolic expression for system dynamics (drift for SDEs). For pure diffusion systems with zero drift, see the detailed Pure Diffusion Systems callout.
```
self._f_sym = sp.Matrix([dC_A_dt, dC_B_dt, dT_dt])
# For pure diffusion: self._f_sym = sp.zeros(n_states, 1)
```
self.diffusion_expr (stochastic systems only): Diffusion matrix for SDE noise structure
```
self.diffusion_expr = sp.Matrix([[sigma_A, 0, 0], ...])
```

Optional Elements with Sensible Defaults

These fields have automatic defaults if not specified:

`self.output_vars` (default: `[]`)

List of computed output symbols (not states). Leave empty if you only observe states.

self.output_vars = []  # Default - no computed outputs

`self._h_sym` (default: full state vector)

Output map $y = h(x)$. If not specified, defaults to observing all states.

# Explicitly set (same as default for our reactor):
self._h_sym = sp.Matrix([C_A, C_B, T])

# If omitted, framework uses: self._h_sym = sp.Matrix(self.state_vars)

`self.order` (default: `1`)

System order defines how you represent system dynamics. Nearly all systems use order=1 (first-order state-space form). See the detailed explanation below for when higher-order forms are useful.

self.order = 1  # Default - first-order state-space form

`self.sde_type` (default: `"ito"`)

For stochastic systems, specifies SDE interpretation.

self.sde_type = "ito"  # Default
# Alternative: self.sde_type = "stratonovich"

Optional Methods

`setup_equilibria()`

If defined, automatically called after system initialization to add equilibria. Useful when equilibria depend on system parameters.

def define_system(self, T_amb_val=300.0, ...):
    # Store parameter value for use in setup_equilibria
    self.T_amb_val = T_amb_val
    
    T_amb = sp.symbols('T_amb', positive=True)
    self.parameters = {T_amb: T_amb_val, ...}
    # ... rest of define_system ...

def setup_equilibria(self):
    # Add parameter-dependent equilibria using stored values
    T_amb = self.T_amb_val
    self.add_equilibrium("complete", x_eq=np.array([0.0, 0.0, T_amb]), ...)

If not defined, only the origin equilibrium is added by default.

Custom Methods and Fields

You can add any additional methods or fields that don’t conflict with base class names:

def define_system(self, C_A0=1.0, T0=350.0, volume=10.0, ...):
    # ... standard setup ...
    
    # Store values needed in custom methods
    self.C_A0 = C_A0  # Initial condition
    self.T0 = T0      # Initial temperature
    self.reactor_volume = volume  # Custom application data
    self.reaction_pathway = "A->B->C"  # Documentation
    
def compute_conversion_rate(self, C_A_current):
    """Calculate conversion from initial concentration."""
    return (self.C_A0 - C_A_current) / self.C_A0

def compute_reactor_capacity(self):
    """Calculate moles of A that can be processed."""
    return self.C_A0 * self.reactor_volume

What NOT to store:

# Don't store state/control symbols as attributes if there are no plans to change parameters post-initialization
# otherwise, can be used to retrieve and re-define paramters in parameter dictionary
# self.C_A_sym = sp.symbols('C_A')  # Redundant!

# Don't store temporary calculations
# self.r1 = k1 * C_A * sp.exp(-E1/T)  # Use local variable

# Don't store parameter values not used elsewhere
# If k1 is only in _f_sym, just put it in self.parameters

For our batch reactor, we explicitly set most optional fields for documentation purposes, but only self.order and self.sde_type are redundant (they match defaults).

Understanding System Order in Detail

System Order: First-Order vs Higher-Order Forms

The self.order field specifies how you represent the system dynamics. There are two mathematically equivalent ways to define systems:

First-Order State-Space Form (order=1, default)

The system is written in first-order form where self._f_sym returns all time derivatives.

Example: Second-order pendulum with states $x = [\theta, \dot{\theta}]$

self.state_vars = [theta, theta_dot]  # Both position and velocity
self._f_sym = sp.Matrix([
    theta_dot,                          # d(theta)/dt
    -(g/L)*sp.sin(theta) - (b/m)*theta_dot + u/m  # d(theta_dot)/dt
])
self.order = 1  # Returns ALL derivatives [d(theta)/dt, d(theta_dot)/dt]

This is the standard state-space form used in control theory. For an $n$-th order physical system with generalized coordinate $q$, the state is:

\[x = [q, \dot{q}, \ddot{q}, \ldots, q^{(n-1)}]\]

and self._f_sym returns the full derivative vector:

\[\frac{dx}{dt} = [\dot{q}, \ddot{q}, \ldots, q^{(n)}]\]

Higher-Order Form (order=n)

The system is written in higher-order form where self._f_sym returns only the highest derivative.

Example: Same pendulum with order=2

self.state_vars = [theta, theta_dot]  # Still both variables
self._f_sym = sp.Matrix([
    -(g/L)*sp.sin(theta) - (b/m)*theta_dot + u/m  # ONLY d²(theta)/dt²
])
self.order = 2  # Returns ONLY highest derivative d²(theta)/dt²

The state vector is the same $x = [\theta, \dot{\theta}]$, but self._f_sym only returns $\ddot{\theta}$. The framework automatically constructs the full first-order representation during linearization by adding the kinematic relationships:

\[\frac{d\theta}{dt} = \dot{\theta}, \quad \frac{d\dot{\theta}}{dt} = \ddot{\theta}\]

When to Use Each Form

First-order (order=1):
- Use this for most systems - it’s explicit and standard
- Required when dynamics aren’t naturally hierarchical
- Example: Our batch reactor has independent first derivatives for $C_A$, $C_B$, and $T$
Higher-order (order=n):
- Natural for mechanical systems (position → velocity → acceleration)
- More compact when you have $\ddot{q} = f(q, \dot{q}, u)$ form
- Example: Robotics, where you directly compute joint accelerations

Requirements

For order=n > 1: The number of states must be divisible by n (i.e., nx % n == 0)
State variables should be ordered as $[q_1, \ldots, q_m, \dot{q}_1, \ldots, \dot{q}_m, \ldots]$ for proper automatic construction
The framework handles all conversions transparently during linearization

This Tutorial

Our batch reactor uses order=1 (default) because the three state variables $(C_A, C_B, T)$ have independent first-order dynamics—there’s no natural hierarchy to exploit.

Equilibrium Setup Example

Our reactor’s setup_equilibria() method demonstrates how to define parameter-dependent equilibria:

def setup_equilibria(self):
    # Use stored parameter values directly
    # (Don't try to look up via sp.symbols() - that creates NEW symbol objects!)
    T_amb = self.T_amb_val

    # Complete conversion equilibrium (deterministic part)
    self.add_equilibrium(
        "complete",
        x_eq=np.array([0.0, 0.0, T_amb]),
        u_eq=np.array([0.0]),
        verify=True,
        stability="stable",
        notes="Equilibrium of deterministic part (drift). Stochastic trajectories "
        "fluctuate around this point with variance growing over time.",
    )

    # Initial condition (if provided)
    if self.C_A0 is not None and self.T0 is not None:
        alpha = self.alpha_val
        Q_init = alpha * (self.T0 - T_amb)

        self.add_equilibrium(
            "initial",
            x_eq=np.array([self.C_A0, 0.0, self.T0]),
            u_eq=np.array([Q_init]),
            verify=False,
            stability="unstable",
            notes="Initial state setpoint (deterministic part). Stochastic trajectories "
            "will fluctuate with std ~ [σ_A·√t, σ_B·√t, σ_T·√t].",
        )
        self.set_default_equilibrium("initial")
    else:
        self.set_default_equilibrium("complete")

Accessing Parameters in setup_equilibria()

Common mistake: Trying to access parameters via self.parameters[sp.symbols("param_name")]

# WRONG: Creates a NEW symbol, won't match dictionary key
T_amb = self.parameters[sp.symbols("T_amb")]  # KeyError!

Why this fails: Each call to sp.symbols() creates a new Python object with unique identity. The dictionary keys in self.parameters are the original symbol objects from define_system(), not new ones.

Solution 1: Store parameter values as attributes (Recommended for beginners):

def define_system(self, T_amb_val=300.0, ...):
    # Store the numerical value
    self.T_amb_val = T_amb_val
    
    # Create symbol and use in parameters dict
    T_amb = sp.symbols('T_amb', positive=True)
    self.parameters = {T_amb: T_amb_val, ...}

def setup_equilibria(self):
    # Use the stored value directly
    T_amb = self.T_amb_val  # Works!

Solution 2: Store symbol references (Recommended for production):

def define_system(self, T_amb_val=300.0, ...):
    # Create symbol and store reference
    self.T_amb_sym = sp.symbols('T_amb', positive=True)
    self.parameters = {self.T_amb_sym: T_amb_val, ...}

def setup_equilibria(self):
    # Use the stored symbol to look up value
    T_amb = self.parameters[self.T_amb_sym]  # Works!

Our tutorial uses Solution 1 for simplicity and clarity.

Simulation and Visualization

Now that we’ve defined our stochastic batch reactor, let’s see it in action. This section demonstrates the complete workflow: inspecting the system, running simulations, and visualizing results.

Integration Methods: integrate() vs simulate()

ControlDESymulation provides two methods for time evolution:

`integrate()` - Low-Level Interface

For stochastic systems, always use this method.

result = system.integrate(
    x0=x0,
    u=controller,
    t_span=(0, 10),
    method='euler_maruyama',  # Stochastic integrator
    dt=0.01,                  # Required for fixed-step methods
    seed=42,                  # Reproducibility
    n_paths=100               # Monte Carlo ensemble
)

Stochastic-specific features:

seed: Random number generator seed for reproducible noise
n_paths: Built-in Monte Carlo simulation (returns shape (n_paths, T, nx))
Stochastic methods: 'euler_maruyama', 'milstein', 'heun', etc.
dt: Time step (required for most SDE methods)
t_eval: Optional array of specific output times

Fixed-Step vs Adaptive Methods:

Most SDE integrators are fixed-step (Euler-Maruyama, Milstein, Heun, etc.):

result = system.integrate(
    x0=x0, t_span=(0, 10),
    method='euler_maruyama',
    dt=0.01  # Required - sets uniform time grid
)
# result['t'] has uniform spacing: [0.00, 0.01, 0.02, ..., 10.00]

Some advanced methods are adaptive (use rtol/atol instead of dt):

result = system.integrate(
    x0=x0, t_span=(0, 10),
    method='adaptive_sde',  # Hypothetical adaptive method
    rtol=1e-6, atol=1e-8
)
# result['t'] may have non-uniform spacing

Visualization tip: For uniform time grids (most cases), plotting is straightforward. For adaptive grids, use t_eval to force specific output times:

result = system.integrate(
    x0=x0, t_span=(0, 10),
    method='adaptive_sde',
    t_eval=np.linspace(0, 10, 1000),  # Force uniform grid
    rtol=1e-6
)

Returns: SDEIntegrationResult with fields:

t: Time points (T,) - uniform for fixed-step, may vary for adaptive
x: States (T, nx) or (n_paths, T, nx) for Monte Carlo
success, nfev, diffusion_evals: Solver diagnostics

`simulate()` - High-Level Interface

For deterministic systems only (not recommended for SDEs).

result = system.simulate(
    x0=x0,
    controller=lambda x, t: -K @ x,
    t_span=(0, 10),
    dt=0.01
)

Key differences:

Always returns regular time grid (uniform spacing dt)
Cleaner API for feedback controllers
No seed or n_paths support
Uses deterministic integrators (RK45, RK4, Euler)
Returns SimulationResult in time-major format (T, nx)

When to use each:

Stochastic systems (SDEs) → integrate() with stochastic methods
Deterministic systems (ODEs) → Either method; simulate() is simpler
Monte Carlo simulations → integrate() with n_paths > 1

For our stochastic batch reactor, we’ll use integrate() exclusively.

Inspecting the System

Before running simulations, let’s examine what we’ve created:

Code

# Display symbolic equations
reactor = ContinuousStochasticBatchReactor(
    C_A0=1.0,    # Initial concentration
    T0=350.0,    # Initial temperature
    sigma_T=2.0  # Temperature noise intensity
)

# Print the mathematical structure
reactor.print_equations(simplify=True)

======================================================================
ContinuousStochasticBatchReactor (Continuous-Time)
======================================================================
State Variables: [C_A, C_B, T]
Control Variables: [Q]
System Order: 1
Dimensions: nx=3, nu=1, ny=3

Dynamics: dx/dt = f(x, u)
  dC_A/dt = -0.5*C_A*exp(-1000.0/T)
  dC_B/dt = 0.5*C_A*exp(-1000.0/T) - 0.3*C_B*exp(-1500.0/T)
  dT/dt = Q - 0.1*T + 30.0

Output: y = h(x)
  y[0] = C_A
  y[1] = C_B
  y[2] = T
======================================================================

This shows:

State and control variables
System dimensions (nx=3 states, nu=1 control)
Symbolic drift dynamics with parameters substituted
Output map (if defined)

We can also inspect system properties programmatically:

Code

print(f"Number of states: {reactor.nx}")
print(f"Number of controls: {reactor.nu}")
print(f"State variables: {reactor.state_vars}")
print(f"Parameters: {list(reactor.parameters.keys())}")

Number of states: 3
Number of controls: 1
State variables: [C_A, C_B, T]
Parameters: [k1, k2, E1, E2, alpha, T_amb, sigma_A, sigma_B, sigma_T]

Choosing an Integration Method

For this tutorial, we use Euler-Maruyama (method='euler_maruyama'), the simplest fixed-step SDE method. It’s a good starting point because:

Fixed time step: Requires dt parameter, produces uniform time grid
Easy visualization: Uniform spacing makes plotting straightforward
Stable for additive noise: Our reactor has additive noise (constant diffusion)
Fast: Minimal computation per step

Other common methods:

milstein: Higher accuracy for multiplicative noise, still fixed-step
heun: Good balance of accuracy and stability, fixed-step
Advanced methods: May use adaptive stepping (see note below)

Handling Adaptive SDE Methods

Some advanced SDE solvers use adaptive time stepping for efficiency. If you use an adaptive method and need uniform output times (e.g., for certain visualizations or statistical analysis), use the t_eval parameter:

# Adaptive method with non-uniform internal steps
result = system.integrate(
    x0=x0,
    t_span=(0, 10),
    method='adaptive_sde_method',  # Hypothetical adaptive method
    rtol=1e-6, atol=1e-8,
    t_eval=np.linspace(0, 10, 1000)  # Force uniform output grid
)

The solver will use adaptive steps internally for efficiency but interpolate to provide output at the exact times specified in t_eval. This gives you:

Efficiency of adaptive stepping during integration
Uniform time grid for downstream analysis and visualization
Control over output density

For this tutorial: We use fixed-step Euler-Maruyama, so dt controls both internal steps and output times. No t_eval needed.

For production code with multiplicative noise, consider higher-order methods. For learning and visualization, Euler-Maruyama with dt=0.01 to dt=0.1 works well.

Equilibrium Points

Our setup_equilibria() method created two equilibrium points. Let’s see them:

Code

# List all equilibria
equilibrium_names = reactor.list_equilibria()
print(f"Available equilibria: {equilibrium_names}")

# Get the "complete conversion" equilibrium
x_eq, u_eq = reactor.get_equilibrium('complete')
print(f"\nComplete conversion equilibrium:")
print(f"  C_A = {x_eq[0]:.3f} mol/L")
print(f"  C_B = {x_eq[1]:.3f} mol/L")
print(f"  T = {x_eq[2]:.1f} K")
print(f"  Q = {u_eq[0]:.3f} W")

# Verify it's actually an equilibrium (drift should be zero)
drift = reactor(x_eq, u_eq, t=0)
print(f"\nDrift at equilibrium: {drift}")
print(f"Max drift magnitude: {np.abs(drift).max():.2e}")

Available equilibria: ['origin', 'complete', 'initial']

Complete conversion equilibrium:
  C_A = 0.000 mol/L
  C_B = 0.000 mol/L
  T = 300.0 K
  Q = 0.000 W

Drift at equilibrium: [-0.  0.  0.]
Max drift magnitude: 0.00e+00

Key insight: For stochastic systems, equilibria are fixed points of the drift only. The diffusion causes trajectories to fluctuate around these points with variance growing as $\sqrt{t}$.

Summary: Building Your System

To define a symbolic system in ControlDESymulation:

Minimum viable system:

def define_system(self):
    # 1. Define symbolic variables
    x = sp.symbols('x')
    u = sp.symbols('u')
    
    # 2. Set required fields
    self.state_vars = [x]
    self.control_vars = [u]  # or [] for autonomous
    self.parameters = {}     # or {param: value}
    self._f_sym = sp.Matrix([...])  # Your dynamics
    
    # 3. For stochastic systems, add:
    # self.diffusion_expr = sp.Matrix([...])

Autonomous system (no control inputs):

def define_system(self, damping=0.1, stiffness=1.0):
    # 1. Define symbolic variables (no control symbols needed)
    x, v = sp.symbols('x v')
    b, k = sp.symbols('b k', positive=True)
    
    # 2. Set required fields
    self.state_vars = [x, v]
    self.control_vars = []  # Empty list for autonomous
    self.parameters = {b: damping, k: stiffness}
    
    # 3. Dynamics depending only on states
    self._f_sym = sp.Matrix([
        v,                  # dx/dt = v
        -k*x - b*v         # dv/dt = -kx - bv (damped oscillator)
    ])

Pure diffusion system (stochastic, zero drift):

def define_system(self, sigma_x=1.0, sigma_y=1.0, sigma_z=1.0):
    # 1. Define symbolic variables
    x, y, z = sp.symbols('x y z', real=True)
    sx, sy, sz = sp.symbols('sigma_x sigma_y sigma_z', positive=True)
    
    # 2. Set required fields
    self.state_vars = [x, y, z]
    self.control_vars = []  # Typically autonomous
    self.parameters = {sx: sigma_x, sy: sigma_y, sz: sigma_z}
    
    # 3. ZERO DRIFT: Pure diffusion has no deterministic dynamics
    self._f_sym = sp.zeros(3, 1)  # Must be zero vector
    
    # 4. DIFFUSION: Stochastic dynamics only
    self.diffusion_expr = sp.Matrix([
        [sx, 0, 0],     # Independent noise in x
        [0, sy, 0],     # Independent noise in y
        [0, 0, sz]      # Independent noise in z
    ])
    
    self.sde_type = "ito"

Well-documented system (recommended):

def define_system(self, param1=1.0, param2=2.0):
    # 1. Define all symbolic variables (kept as LOCAL variables)
    x, y = sp.symbols('x y')        # State symbols - local only
    u = sp.symbols('u')             # Control symbol - local only
    p1, p2 = sp.symbols('p1 p2', positive=True)  # Parameter symbols - local only
    
    # 2. Set required fields
    self.state_vars = [x, y]        # Reference to local symbols
    self.control_vars = [u]         # Reference to local symbols
    self.parameters = {p1: param1, p2: param2}  # Symbol→value mapping
    self._f_sym = sp.Matrix([...])
    
    # 3. Explicitly set optional fields for clarity
    self.output_vars = []           # No computed outputs
    self._h_sym = sp.Matrix([x, y]) # Observe full state
    self.order = 1                  # First-order form
    
    # 4. For stochastic systems:
    # self.diffusion_expr = sp.Matrix([...])
    # self.sde_type = "ito"  # or "stratonovich"
    
    # 5. Store parameter VALUES as attributes (if needed in setup_equilibria)
    # DON'T store symbols unless parameters need to change post-initialization
    # symbols are already in self.state_vars, self.parameters
    # self.param1_val = param1  # Store VALUE for later use
    # self.param2_val = param2  # Store VALUE for later use

def setup_equilibria(self):
    # Optional: Add parameter-dependent equilibria
    # Access stored values: x_eq_value = self.param1_val
    self.add_equilibrium("my_equilibrium", x_eq=..., u_eq=...)

Which template should you use?

Minimum viable: Quick prototyping, testing, simple demos
Autonomous system: No control inputs (free oscillators, population models, unforced reactions)
Pure diffusion: Stochastic processes with zero drift (Brownian motion, pure noise)
Well-documented: Production code, complex systems, collaborative projects (recommended)

The reactor example in this tutorial follows the well-documented approach with both drift and diffusion.

What You’ve Learned

This tutorial covered the anatomy of a symbolic system through a stochastic batch reactor example:

System Definition

How to subclass ContinuousStochasticSystem
Required fields: state_vars, control_vars, parameters, _f_sym, diffusion_expr
Optional fields: output_vars, _h_sym, order, sde_type
Special cases: Autonomous systems and pure diffusion systems
Parameter-dependent equilibria via setup_equilibria()
What to store as attributes vs local variables:
- Symbols (state, control, parameters) → Keep local
- Parameter values (used in setup_equilibria) → Store as self.param_val
- Initial conditions and custom data → Store as attributes

System Introspection

print_equations(): Display symbolic dynamics
list_equilibria(): See all equilibrium points
get_equilibrium(): Access equilibrium states and controls
Programmatic access: nx, nu, state_vars, parameters

Critical Distinctions

Drift vs Diffusion: Deterministic dynamics vs stochastic forcing
Itô vs Stratonovich: Different stochastic calculus interpretations
Equilibria in SDEs: Fixed points of drift only (stochastic fluctuations persist)
integrate() vs simulate(): Use integrate() for stochastic systems
Fixed-step vs adaptive: Most SDE methods use fixed dt (uniform time grid), some advanced methods adapt step size (use t_eval for uniform output)
Parameter access: Store parameter values as attributes (self.param_val) for use in setup_equilibria(), don’t try to look up via self.parameters[sp.symbols("param")] (creates new symbol, causes KeyError). Store symbol references as attributes only if parameter values need to change post-initialization.

Next Steps

Tutorial 2: Simulation and Visualization

Single trajectory simulation with integrate()
Monte Carlo ensembles with n_paths
Time series visualization with reactor.plotter
Phase portraits and state-space analysis
Ensemble statistics and confidence bands

--- title: "Anatomy of a Symbolic System" subtitle: "Stochastic Batch Reactors" author: "Gil Benezer" date: today format: html: toc: true toc-depth: 5 code-fold: show code-tools: true theme: cosmo execute: eval: true cache: true warning: false --- This tutorial introduces the anatomy of a symbolic control system using a **stochastic batch reactor** as a worked example. A batch reactor is a closed chemical system where reactants $A$ and $B$ undergo sequential reactions governed by Arrhenius kinetics, with temperature controlled via external heating $Q$. ## System Overview The system demonstrates key `cdesym` concepts: - **State variables**: Concentrations $(C_A, C_B)$ and temperature $(T)$ - **Control input**: Heat input $(Q)$ - **Drift dynamics**: Deterministic reaction rates and heat transfer - **Diffusion dynamics**: Additive noise on each state, representing measurement or process uncertainty The reaction sequence $A \to B \to C$ follows first-order Arrhenius kinetics: $$ \begin{aligned} r_1 &= k_1 \cdot C_A \cdot \exp(\frac{-E_1}{T}) \\ r_2 &= k_2 \cdot C_B \cdot \exp(\frac{-E_2}{T}) \end{aligned} $$ The deterministic dynamics are: $$ \begin{aligned} \frac{dC_A}{dt} &= -r_1 \\ \frac{dC_B}{dt} &= r_1 - r_2 \\ \frac{dT}{dt} &= Q - \alpha(T - T_{\text{amb}}) \end{aligned} $$ ::: {.callout-note collapse="true"} ## Itô vs Stratonovich SDEs The stochastic formulation in this example uses an **Itô SDE** of the form: $$dX_t = f(X_t, u_t)\,dt + g(X_t)\,dW_t$$ where $W_t$ is a Wiener process (Brownian motion). The key distinction from **Stratonovich SDEs** lies in how the stochastic integral is defined: | Property | Itô | Stratonovich | |----------|-----|--------------| | **Evaluation point** | Left endpoint of interval | Midpoint of interval | | **Martingale property** | Yes (integral is a martingale) | No | | **Chain rule** | Requires **Itô's lemma** with correction term | Standard chain rule applies | | **Typical use** | Finance, control theory, probability | Physics, systems with colored noise limits | **Itô's lemma** for a function $Y_t = h(X_t)$ includes an extra term: $$dY_t = h'(X_t)\,dX_t + \frac{1}{2}h''(X_t)\,g(X_t)^2\,dt$$ The second term (absent in ordinary calculus) arises because Brownian paths have infinite quadratic variation. **When the distinction matters**: For **state-dependent (multiplicative) noise** like $g(X_t) = \sigma X_t$, the two interpretations give different results. For **additive noise** where $g$ is constant, they coincide. This tutorial uses additive noise, so the choice is moot here—but understanding the difference is essential when modeling state-dependent noise. ::: ## Noise Structure The stochastic batch reactor extends the deterministic model with additive process noise: $$ \begin{aligned} dC_A &= (-r_1)\,dt + \sigma_A\,dW_A \\ dC_B &= (r_1 - r_2)\,dt + \sigma_B\,dW_B \\ dT &= (Q - \alpha(T - T_{\text{amb}}))\,dt + \sigma_T\,dW_T \end{aligned} $$ Three independent Wiener processes model distinct physical noise sources: feed variability ($\sigma_A$, $\sigma_B$) and heat transfer fluctuations ($\sigma_T$). Trajectories fluctuate around the deterministic equilibrium with standard deviation growing as $\sqrt{t}$. ## What This Tutorial Covers We define two equilibria: 1. **Complete conversion** (stable): All reactants consumed, temperature at ambient 2. **Initial setpoint** (unstable): Starting conditions for batch operation The tutorial walks through defining `_f_sym` (drift), `diffusion_expr` (noise structure), how to define optional `_h_sym` (output map) and output variables if needed, then covers equilibrium setup and trajectory visualization. ## System Definition The full system definition is as follows: ```{python} from typing import Optional import numpy as np import sympy as sp from cdesym import ContinuousStochasticSystem class ContinuousStochasticBatchReactor(ContinuousStochasticSystem): def define_system( self, k1_val: float = 0.5, k2_val: float = 0.3, E1_val: float = 1000.0, E2_val: float = 1500.0, alpha_val: float = 0.1, T_amb_val: float = 300.0, sigma_A: float = 0.01, sigma_B: float = 0.01, sigma_T: float = 1.0, C_A0: Optional[float] = None, T0: Optional[float] = None, ): # Store initial conditions and parameter values for later use # (See "What to Store as Attributes" callout below for detailed guidelines) self.C_A0 = C_A0 self.T0 = T0 self.alpha_val = alpha_val self.T_amb_val = T_amb_val # State and control variables (local - not stored as attributes) C_A, C_B, T = sp.symbols("C_A C_B T", real=True, positive=True) Q = sp.symbols("Q", real=True) # Deterministic Parameters (kinetics and heat transfer) # These are local variables - used as keys in self.parameters dict k1, k2, E1, E2, alpha, T_amb = sp.symbols( "k1 k2 E1 E2 alpha T_amb", real=True, positive=True, ) # Stochastic Parameters (noise intensities) sigma_A_sym = sp.symbols("sigma_A", real=True, positive=True) sigma_B_sym = sp.symbols("sigma_B", real=True, positive=True) sigma_T_sym = sp.symbols("sigma_T", real=True, positive=True) # Parameters dictionary: Symbol → Value mapping # Keys are the symbol objects created above self.parameters = { k1: k1_val, k2: k2_val, E1: E1_val, E2: E2_val, alpha: alpha_val, T_amb: T_amb_val, sigma_A_sym: sigma_A, sigma_B_sym: sigma_B, sigma_T_sym: sigma_T, } self.state_vars = [C_A, C_B, T] self.control_vars = [Q] self.output_vars = [] self.order = 1 # Reaction rates (Arrhenius kinetics) r1 = k1 * C_A * sp.exp(-E1 / T) r2 = k2 * C_B * sp.exp(-E2 / T) # DRIFT (Deterministic part - same as deterministic reactor) dC_A_dt = -r1 dC_B_dt = r1 - r2 dT_dt = Q - alpha * (T - T_amb) self._f_sym = sp.Matrix([dC_A_dt, dC_B_dt, dT_dt]) # DIFFUSION (Stochastic part - additive noise) # Diagonal matrix: three independent Wiener processes self.diffusion_expr = sp.Matrix( [ [sigma_A_sym, 0, 0], [0, sigma_B_sym, 0], [0, 0, sigma_T_sym], ], ) # Itô SDE interpretation self.sde_type = "ito" # Output: Full state measurement (with potential noise in practice) # technically optional self._h_sym = sp.Matrix([C_A, C_B, T]) def setup_equilibria(self): # Use stored parameter values directly # (Don't try to look up via sp.symbols() - that creates NEW symbol objects!) T_amb = self.T_amb_val # Complete conversion equilibrium (deterministic part) self.add_equilibrium( "complete", x_eq=np.array([0.0, 0.0, T_amb]), u_eq=np.array([0.0]), verify=True, stability="stable", notes="Equilibrium of deterministic part (drift). Stochastic trajectories " "fluctuate around this point with variance growing over time.", ) # Initial condition (if provided) if self.C_A0 is not None and self.T0 is not None: alpha = self.alpha_val Q_init = alpha * (self.T0 - T_amb) self.add_equilibrium( "initial", x_eq=np.array([self.C_A0, 0.0, self.T0]), u_eq=np.array([Q_init]), verify=False, stability="unstable", notes="Initial state setpoint (deterministic part). Stochastic trajectories " "will fluctuate with std ~ [σ_A·√t, σ_B·√t, σ_T·√t].", ) self.set_default_equilibrium("initial") else: self.set_default_equilibrium("complete") ``` The first step in system construction is to subclass the proper symbolic system and create the define_system method signature. ```python class ContinuousStochasticBatchReactor(ContinuousStochasticSystem): def define_system( self, k1_val: float = 0.5, k2_val: float = 0.3, E1_val: float = 1000.0, E2_val: float = 1500.0, alpha_val: float = 0.1, T_amb_val: float = 300.0, sigma_A: float = 0.01, sigma_B: float = 0.01, sigma_T: float = 1.0, C_A0: Optional[float] = None, T0: Optional[float] = None, ): ``` More often than not, the only arguments to this function will be the `self` argument and parameter values. It is not necessary to add type hints or default values, but it is recommended for debug purposes. The next step in system definition is to define the (deterministic) symbolic variables to be used and place them in the proper fields. The main limitation currently is that it is not possible to explicitly define a symbolic variable for time. If it is heavily requested, we can implement it as a future feature. ```python # State and control variables C_A, C_B, T = sp.symbols("C_A C_B T", real=True, positive=True) Q = sp.symbols("Q", real=True) # Deterministic Parameters (kinetics and heat transfer) k1, k2, E1, E2, alpha, T_amb = sp.symbols( "k1 k2 E1 E2 alpha T_amb", real=True, positive=True, ) # Stochastic Parameters (noise intensities) sigma_A_sym = sp.symbols("sigma_A", real=True, positive=True) sigma_B_sym = sp.symbols("sigma_B", real=True, positive=True) sigma_T_sym = sp.symbols("sigma_T", real=True, positive=True) self.parameters = { k1: k1_val, k2: k2_val, E1: E1_val, E2: E2_val, alpha: alpha_val, T_amb: T_amb_val, sigma_A_sym: sigma_A, sigma_B_sym: sigma_B, sigma_T_sym: sigma_T, } self.state_vars = [C_A, C_B, T] self.control_vars = [Q] self.output_vars = [] ``` The necessary symbolic variables include: - State variables - Control variables - Parameter variables ::: {.callout-warning} ## What to Store as Attributes vs Local Variables A common source of confusion: **when should you store something as `self.something` vs keep it as a local variable?** ### SymPy Symbols: Usually Local, Sometimes Stored **State and control symbols → Keep LOCAL (don't store as attributes):** ```python # CORRECT: Local variables only C_A, C_B, T = sp.symbols("C_A C_B T", real=True) Q = sp.symbols("Q", real=True) # WRONG: Don't store these as attributes # self.C_A_sym = sp.symbols("C_A") # Unnecessary! ``` **Why keep local?** These symbols are only used to build `self._f_sym` and `self._h_sym`. Once those are built, the framework tracks symbols via `self.state_vars` and `self.control_vars` lists. Storing individual symbols is redundant. **Parameter symbols → Keep LOCAL (unless using Solution 2 for equilibria):** ```python # Solution 1 (tutorial approach): Local symbols T_amb = sp.symbols('T_amb', positive=True) self.parameters = {T_amb: T_amb_val} # Symbol used as key only # Solution 2 (production code): Store symbol references self.T_amb_sym = sp.symbols('T_amb', positive=True) self.parameters = {self.T_amb_sym: T_amb_val} # Can look up later ``` ### Parameter Values: Store When Needed Later **Store parameter values if used in `setup_equilibria()` or custom methods:** ```python def define_system(self, T_amb_val=300.0, alpha_val=0.1, ...): # Store values needed in setup_equilibria self.T_amb_val = T_amb_val self.alpha_val = alpha_val # Create symbols and populate parameters dict T_amb = sp.symbols('T_amb', positive=True) self.parameters = {T_amb: T_amb_val, ...} ``` **Don't store if only used in symbolic expressions:** ```python def define_system(self, k1_val=0.5, k2_val=0.3, ...): # These are only used in _f_sym, not in setup_equilibria # No need to store as attributes k1 = sp.symbols('k1', positive=True) k2 = sp.symbols('k2', positive=True) self.parameters = {k1: k1_val, k2: k2_val} self._f_sym = sp.Matrix([k1 * C_A * sp.exp(-E1/T)]) # Used here only ``` ### Custom Fields: Always Store **Store data needed for conditional logic or custom methods:** ```python def define_system(self, C_A0=None, T0=None, ...): # Store for use in setup_equilibria self.C_A0 = C_A0 self.T0 = T0 # Store custom application data self.reaction_pathway = "A->B->C" self.reactor_volume = 10.0 # liters ``` ### Decision Tree **Should I store `foo` as `self.foo`?** 1. **Is it a SymPy symbol?** - State/control symbol? → **No** (use `self.state_vars` / `self.control_vars` lists) - Parameter symbol? → **No** (unless using Solution 2 for parameter lookup) 2. **Is it a parameter value?** - Used in `setup_equilibria()`? → **Yes** (`self.param_val`) - Only used in `_f_sym`? → **No** (just put in `self.parameters` dict) 3. **Is it a custom field?** - Used in logic or custom methods? → **Yes** (`self.C_A0`, `self.reactor_volume`) - Temporary calculation? → **No** (keep as local variable) **Rule of thumb:** If you'll need it in `setup_equilibria()` or a custom method, store it as `self.something`. Otherwise, keep it local. ::: ::: {.callout-warning} ## Output Variables and the Output Map There are two related but distinct concepts for system outputs: - **`self.output_vars`**: A list of symbolic variables representing *computed outputs*—quantities derived from state variables but not states themselves. For example, if your states are position and velocity, an output variable might be kinetic energy. Do **not** add state variables to this list. - **`self._h_sym`**: The *output map* $y = h(x)$, a symbolic column vector defining what the system actually outputs. This can include: - State variables (for full or partial state observation) - Output variables (computed quantities) - Any symbolic expression of states **Defaults**: `self.output_vars = []` (no computed outputs) and `self._h_sym = state vector` (full state observability). **Example**: For a system with states $[x, v]$ and a computed output $E = \frac{1}{2}mv^2$: - `self.output_vars = [E]` - `self._h_sym = sp.Matrix([x, E])` outputs position and energy (but not velocity) ::: ::: {.callout-note collapse="true"} ## Autonomous Systems An **autonomous system** has no external control inputs—its dynamics depend only on the current state (and possibly time). Examples include free oscillators, population dynamics, and unforced chemical reactions. ### Required Specification To define an autonomous system, you **must** set `self.control_vars = []` (an empty list). Do not use `None` or omit the field—this will cause errors during system initialization. ```python # CORRECT: Autonomous system self.control_vars = [] # WRONG: Will cause errors self.control_vars = None # self.control_vars = ... (omitting the field) ``` ### Dynamics for Autonomous Systems For autonomous systems, `self._f_sym` should depend only on state variables and parameters, **not control variables**: ```python # Example: Damped harmonic oscillator (autonomous) x, v = sp.symbols('x v', real=True) k, b, m = sp.symbols('k b m', positive=True) self.state_vars = [x, v] self.control_vars = [] # No control self.parameters = {k: 10.0, b: 0.5, m: 1.0} # Dynamics: dx/dt = v, dv/dt = -(k/m)*x - (b/m)*v self._f_sym = sp.Matrix([ v, -(k/m)*x - (b/m)*v ]) ``` ### Calling Conventions When evaluating or integrating autonomous systems: - The `u` parameter in `system(x, u)` should be `None` or omitted - Integration accepts `u=None` or no control specification - The framework automatically handles zero control internally ```python # Evaluating dynamics dxdt = system(x) # Recommended dxdt = system(x, u=None) # Also valid dxdt = system(x, np.array([])) # Works but unnecessary # Integration result = system.integrate(x0=x_init, t_span=(0, 10)) # No need to specify u - it's ignored for autonomous systems ``` ### Converting Between Controlled and Autonomous You can model the same physical system in different ways: **Controlled oscillator** (with external force): ```python self.state_vars = [x, v] self.control_vars = [u] # External force self._f_sym = sp.Matrix([v, -(k/m)*x - (b/m)*v + u/m]) ``` **Free oscillator** (autonomous, special case u=0): ```python self.state_vars = [x, v] self.control_vars = [] # No control self._f_sym = sp.Matrix([v, -(k/m)*x - (b/m)*v]) ``` Both are valid modeling choices. Use controlled form when you plan to apply inputs at some point; use autonomous form for analyzing free dynamics or when control is not part of the problem. Systems with control inputs expect at least a function that output constant zero input. ### Time-Varying vs Autonomous **Important distinction**: - **Autonomous** = no control inputs (`self.control_vars = []`) - **Time-invariant** = dynamics don't explicitly depend on time - **Time-varying** = dynamics explicitly depend on time A system can be autonomous but time-varying (e.g., $\dot{x} = -\sin(t) \cdot x$) or controlled but time-invariant (e.g., $\dot{x} = -x + u$). Most systems in control theory are both controlled and time-invariant. ### Practical Examples **Autonomous systems:** - Lotka-Volterra predator-prey: $\dot{x} = \alpha x - \beta xy$, $\dot{y} = \delta xy - \gamma y$ - Van der Pol oscillator: $\ddot{x} + \mu(x^2-1)\dot{x} + x = 0$ - Free batch reactor (no heating): $\dot{C}_A = -k_1 C_A e^{-E_1/T}$, $\dot{T} = -\alpha(T - T_{amb})$ **Controlled systems:** - Inverted pendulum with torque: $\ddot{\theta} = (g/L)\sin\theta + u/(mL^2)$ - Batch reactor with heating (this tutorial): $\dot{T} = Q - \alpha(T - T_{amb})$ - Linear regulator: $\dot{x} = Ax + Bu$ ::: ::: {.callout-important} ## Parameter Dictionary Keys Must Be Symbolic Variables A common mistake is to set the dictionary keys for self.parameters as strings, but for proper symbolic substitution and code generation, the keys must be SymPy Symbolic variables. ::: The final required step in system definition is to relate symbolic variables for the definition of the system (and possibly observation) dynamics. ```python # Reaction rates (Arrhenius kinetics) r1 = k1 * C_A * sp.exp(-E1 / T) r2 = k2 * C_B * sp.exp(-E2 / T) # DRIFT (Deterministic part - same as deterministic reactor) dC_A_dt = -r1 dC_B_dt = r1 - r2 dT_dt = Q - alpha * (T - T_amb) self._f_sym = sp.Matrix([dC_A_dt, dC_B_dt, dT_dt]) # DIFFUSION (Stochastic part - additive noise) # Diagonal matrix: three independent Wiener processes self.diffusion_expr = sp.Matrix( [ [sigma_A_sym, 0, 0], [0, sigma_B_sym, 0], [0, 0, sigma_T_sym], ], ) ``` ::: {.callout-note collapse="true"} ## Diffusion Matrix Structure The diffusion matrix $G$ has dimensions $(n_{\text{states}} \times n_{\text{Wiener}})$, where each column corresponds to an independent Wiener process and each row to a state variable. The stochastic term in the SDE is $G\,dW_t$ where $dW_t$ is a vector of independent Wiener increments. **Why a matrix, not a vector?** A vector would only allow one noise source affecting all states with fixed relative magnitudes. The matrix form enables: - **Diagonal matrices** (as above): Each state has its own independent noise source. The $i$-th Wiener process affects only the $i$-th state. - **Off-diagonal elements**: Represent *correlated noise*—when a single physical noise source affects multiple states. For example, if temperature fluctuations also induced concentration measurement errors: ```python self.diffusion_expr = sp.Matrix([ [sigma_A_sym, 0, sigma_AT], # C_A affected by W_A and W_T [0, sigma_B_sym, 0], [0, 0, sigma_T_sym], ]) ``` - **Non-square matrices**: You can have fewer Wiener processes than states (correlated noise) or more (redundant noise sources for modeling flexibility). **Dimensions expected**: For $n$ state variables and $m$ independent Wiener processes, `diffusion_expr` should be an $n \times m$ SymPy Matrix. ::: ```python # Output: Full state measurement (with potential noise in practice) # technically optional self._h_sym = sp.Matrix([C_A, C_B, T]) ``` ::: {.callout-note collapse="true"} ## Pure Diffusion Systems A **pure diffusion system** (also called a pure noise process) has no deterministic drift—dynamics are driven entirely by stochastic forcing. Examples include Brownian motion, geometric random walks, and models of measurement noise. ### Mathematical Form For pure diffusion systems, the SDE has zero drift: $$dX_t = 0 \cdot dt + G(X_t)\,dW_t = G(X_t)\,dW_t$$ This means the state evolution is governed entirely by the diffusion term. The deterministic dynamics $f(x,u) = 0$ vanish. ### Required Specification To define a pure diffusion system, you **must** set `self._f_sym` to a zero vector with the same dimension as the number of state variables. Use SymPy's `zeros()` function or explicitly construct the zero matrix: ```python # CORRECT: Pure diffusion system (3 states) self._f_sym = sp.zeros(3, 1) # Recommended approach # ALSO CORRECT: Explicit zero matrix self._f_sym = sp.Matrix([0, 0, 0]) # ALSO CORRECT: Using symbols and expressions that equal zero self._f_sym = sp.Matrix([x - x, y - y, z - z]) # Works but unnecessary # WRONG: Omitting _f_sym entirely (required field) # self._f_sym = None # Will cause errors ``` ### Example: 3D Brownian Motion A particle undergoing pure diffusion in 3D space: ```python class BrownianMotion3D(ContinuousStochasticSystem): def define_system(self, sigma_x=1.0, sigma_y=1.0, sigma_z=1.0): # Position coordinates x, y, z = sp.symbols('x y z', real=True) # Diffusion coefficients sx, sy, sz = sp.symbols('sigma_x sigma_y sigma_z', positive=True) self.state_vars = [x, y, z] self.control_vars = [] # Autonomous (no control) self.parameters = {sx: sigma_x, sy: sigma_y, sz: sigma_z} # ZERO DRIFT: No deterministic motion self._f_sym = sp.zeros(3, 1) # DIFFUSION: Independent noise in each direction self.diffusion_expr = sp.Matrix([ [sx, 0, 0], [0, sy, 0], [0, 0, sz] ]) self.sde_type = "ito" ``` ### When to Use Pure Diffusion **Use pure diffusion when modeling:** - **Brownian motion**: Pure random walks in continuous time - **Measurement noise**: Sensors with zero mean noise and no drift - **Diffusion-only processes**: Heat diffusion at steady-state, particle diffusion - **Null models**: Testing whether observed dynamics are purely random - **Noise characterization**: Estimating noise covariance without modeling dynamics **Don't use pure diffusion when:** - System has deterministic dynamics (even if noisy) - Modeling physical systems with forces/flows (use drift + diffusion) - Equilibrium points matter (pure diffusion has no equilibria except possibly at boundaries) ### Physical Interpretation For additive noise (constant $G$): $$dX_t = G\,dW_t$$ The solution is simply $X_t = X_0 + G\,W_t$, where $W_t$ is Brownian motion. The state performs a random walk around its initial condition with variance growing linearly: $\text{Var}(X_t) = GG^T \cdot t$. For multiplicative noise (state-dependent $G(X)$): $$dX_t = G(X_t)\,dW_t$$ The dynamics are more complex, and the SDE interpretation (Itô vs Stratonovich) becomes critical. ### Integration Considerations Pure diffusion systems require: - **Stochastic integrators**: Methods like Euler-Maruyama, Milstein, or higher-order SDE schemes - **Fine time steps**: To resolve the noise structure accurately - **Multiple realizations**: Statistical properties emerge from ensemble averaging ```python # Example: Simulating 3D Brownian motion brownian = BrownianMotion3D(sigma_x=1.0, sigma_y=1.0, sigma_z=1.0) # Initial position at origin x0 = np.array([0.0, 0.0, 0.0]) # Integrate using Euler-Maruyama result = brownian.integrate( x0=x0, t_span=(0, 10), method="euler_maruyama", dt=0.01, # Fine time step for noise resolution seed=42 # Reproducible noise ) # For statistics, run multiple trajectories trajectories = [ brownian.integrate(x0=x0, t_span=(0, 10), dt=0.01, seed=i) for i in range(1000) ] ``` ### Drift + Diffusion vs Pure Diffusion Most physical systems have **both** drift and diffusion: **Drift + Diffusion** (this tutorial's reactor): ```python # Deterministic part (drift) self._f_sym = sp.Matrix([dC_A_dt, dC_B_dt, dT_dt]) # Stochastic part (diffusion) self.diffusion_expr = sp.Matrix([[σ_A, 0, 0], ...]) ``` **Pure Diffusion** (rare in practice): ```python # No drift self._f_sym = sp.zeros(3, 1) # Only diffusion self.diffusion_expr = sp.Matrix([[σ_A, 0, 0], ...]) ``` The batch reactor in this tutorial has **drift + diffusion** because chemical reactions (drift) occur alongside process noise (diffusion). ### Practical Applications **Pure Brownian Motion:** - Modeling colloidal particle motion in fluids - Financial models with zero expected return - Null hypothesis testing in time series **State-Dependent Diffusion:** - Volatility modeling in finance (stochastic volatility models) - Birth-death processes at steady-state - Fluctuation-dissipation in statistical mechanics **General Recommendation**: If you're unsure whether your system should be pure diffusion, it probably shouldn't be. Most physical systems have deterministic dynamics (drift) plus noise (diffusion). ::: ## Required vs Optional System Definition Elements Understanding what's required versus optional helps you define systems efficiently while maintaining clarity. Here's the complete breakdown: ::: {.callout-important} ## Required Elements Every `define_system()` method **must** define: 1. **`self.state_vars`**: List of state variable symbols ```python self.state_vars = [C_A, C_B, T] ``` 2. **`self.control_vars`**: List of control variable symbols (empty list `[]` for autonomous systems—see the detailed Autonomous Systems callout) ```python self.control_vars = [Q] # or [] for autonomous ``` 3. **`self.parameters`**: Dictionary mapping symbolic parameters to numerical values ```python self.parameters = {k1: 0.5, k2: 0.3, ...} ``` 4. **`self._f_sym`**: Symbolic expression for system dynamics (drift for SDEs). For pure diffusion systems with zero drift, see the detailed Pure Diffusion Systems callout. ```python self._f_sym = sp.Matrix([dC_A_dt, dC_B_dt, dT_dt]) # For pure diffusion: self._f_sym = sp.zeros(n_states, 1) ``` 5. **`self.diffusion_expr`** (stochastic systems only): Diffusion matrix for SDE noise structure ```python self.diffusion_expr = sp.Matrix([[sigma_A, 0, 0], ...]) ``` ::: ::: {.callout-tip collapse="true"} ## Optional Elements with Sensible Defaults These fields have **automatic defaults** if not specified: ### `self.output_vars` (default: `[]`) List of computed output symbols (not states). Leave empty if you only observe states. ```python self.output_vars = [] # Default - no computed outputs ``` ### `self._h_sym` (default: full state vector) Output map $y = h(x)$. If not specified, defaults to observing all states. ```python # Explicitly set (same as default for our reactor): self._h_sym = sp.Matrix([C_A, C_B, T]) # If omitted, framework uses: self._h_sym = sp.Matrix(self.state_vars) ``` ### `self.order` (default: `1`) System order defines how you represent system dynamics. Nearly all systems use `order=1` (first-order state-space form). See the detailed explanation below for when higher-order forms are useful. ```python self.order = 1 # Default - first-order state-space form ``` ### `self.sde_type` (default: `"ito"`) For stochastic systems, specifies SDE interpretation. ```python self.sde_type = "ito" # Default # Alternative: self.sde_type = "stratonovich" ``` ::: ::: {.callout-tip collapse="true"} ## Optional Methods ### `setup_equilibria()` If defined, automatically called after system initialization to add equilibria. Useful when equilibria depend on system parameters. ```python def define_system(self, T_amb_val=300.0, ...): # Store parameter value for use in setup_equilibria self.T_amb_val = T_amb_val T_amb = sp.symbols('T_amb', positive=True) self.parameters = {T_amb: T_amb_val, ...} # ... rest of define_system ... def setup_equilibria(self): # Add parameter-dependent equilibria using stored values T_amb = self.T_amb_val self.add_equilibrium("complete", x_eq=np.array([0.0, 0.0, T_amb]), ...) ``` If not defined, only the origin equilibrium is added by default. ### Custom Methods and Fields You can add any additional methods or fields that don't conflict with base class names: ```python def define_system(self, C_A0=1.0, T0=350.0, volume=10.0, ...): # ... standard setup ... # Store values needed in custom methods self.C_A0 = C_A0 # Initial condition self.T0 = T0 # Initial temperature self.reactor_volume = volume # Custom application data self.reaction_pathway = "A->B->C" # Documentation def compute_conversion_rate(self, C_A_current): """Calculate conversion from initial concentration.""" return (self.C_A0 - C_A_current) / self.C_A0 def compute_reactor_capacity(self): """Calculate moles of A that can be processed.""" return self.C_A0 * self.reactor_volume ``` **What NOT to store:** ```python # Don't store state/control symbols as attributes if there are no plans to change parameters post-initialization # otherwise, can be used to retrieve and re-define paramters in parameter dictionary # self.C_A_sym = sp.symbols('C_A') # Redundant! # Don't store temporary calculations # self.r1 = k1 * C_A * sp.exp(-E1/T) # Use local variable # Don't store parameter values not used elsewhere # If k1 is only in _f_sym, just put it in self.parameters ``` ::: For our batch reactor, we **explicitly set** most optional fields for documentation purposes, but only `self.order` and `self.sde_type` are redundant (they match defaults). ### Understanding System Order in Detail ::: {.callout-note collapse="true"} ## System Order: First-Order vs Higher-Order Forms The `self.order` field specifies how you represent the system dynamics. There are two mathematically equivalent ways to define systems: ### First-Order State-Space Form (order=1, **default**) The system is written in first-order form where `self._f_sym` returns **all** time derivatives. **Example: Second-order pendulum** with states $x = [\theta, \dot{\theta}]$ ```python self.state_vars = [theta, theta_dot] # Both position and velocity self._f_sym = sp.Matrix([ theta_dot, # d(theta)/dt -(g/L)*sp.sin(theta) - (b/m)*theta_dot + u/m # d(theta_dot)/dt ]) self.order = 1 # Returns ALL derivatives [d(theta)/dt, d(theta_dot)/dt] ``` This is the **standard state-space form** used in control theory. For an $n$-th order physical system with generalized coordinate $q$, the state is: $$x = [q, \dot{q}, \ddot{q}, \ldots, q^{(n-1)}]$$ and `self._f_sym` returns the full derivative vector: $$\frac{dx}{dt} = [\dot{q}, \ddot{q}, \ldots, q^{(n)}]$$ ### Higher-Order Form (order=n) The system is written in higher-order form where `self._f_sym` returns **only the highest** derivative. **Example: Same pendulum** with `order=2` ```python self.state_vars = [theta, theta_dot] # Still both variables self._f_sym = sp.Matrix([ -(g/L)*sp.sin(theta) - (b/m)*theta_dot + u/m # ONLY d²(theta)/dt² ]) self.order = 2 # Returns ONLY highest derivative d²(theta)/dt² ``` The state vector is the same $x = [\theta, \dot{\theta}]$, but `self._f_sym` only returns $\ddot{\theta}$. The framework **automatically constructs** the full first-order representation during linearization by adding the kinematic relationships: $$\frac{d\theta}{dt} = \dot{\theta}, \quad \frac{d\dot{\theta}}{dt} = \ddot{\theta}$$ ### When to Use Each Form - **First-order (order=1)**: - **Use this for most systems** - it's explicit and standard - Required when dynamics aren't naturally hierarchical - Example: Our batch reactor has independent first derivatives for $C_A$, $C_B$, and $T$ - **Higher-order (order=n)**: - Natural for mechanical systems (position → velocity → acceleration) - More compact when you have $\ddot{q} = f(q, \dot{q}, u)$ form - Example: Robotics, where you directly compute joint accelerations ### Requirements - For `order=n > 1`: The number of states must be divisible by `n` (i.e., `nx % n == 0`) - State variables should be ordered as $[q_1, \ldots, q_m, \dot{q}_1, \ldots, \dot{q}_m, \ldots]$ for proper automatic construction - The framework handles all conversions transparently during linearization ### This Tutorial Our batch reactor uses **order=1** (default) because the three state variables $(C_A, C_B, T)$ have independent first-order dynamics—there's no natural hierarchy to exploit. ::: ## Equilibrium Setup Example Our reactor's `setup_equilibria()` method demonstrates how to define parameter-dependent equilibria: ```python def setup_equilibria(self): # Use stored parameter values directly # (Don't try to look up via sp.symbols() - that creates NEW symbol objects!) T_amb = self.T_amb_val # Complete conversion equilibrium (deterministic part) self.add_equilibrium( "complete", x_eq=np.array([0.0, 0.0, T_amb]), u_eq=np.array([0.0]), verify=True, stability="stable", notes="Equilibrium of deterministic part (drift). Stochastic trajectories " "fluctuate around this point with variance growing over time.", ) # Initial condition (if provided) if self.C_A0 is not None and self.T0 is not None: alpha = self.alpha_val Q_init = alpha * (self.T0 - T_amb) self.add_equilibrium( "initial", x_eq=np.array([self.C_A0, 0.0, self.T0]), u_eq=np.array([Q_init]), verify=False, stability="unstable", notes="Initial state setpoint (deterministic part). Stochastic trajectories " "will fluctuate with std ~ [σ_A·√t, σ_B·√t, σ_T·√t].", ) self.set_default_equilibrium("initial") else: self.set_default_equilibrium("complete") ``` ::: {.callout-warning} ## Accessing Parameters in `setup_equilibria()` **Common mistake:** Trying to access parameters via `self.parameters[sp.symbols("param_name")]` ```python # WRONG: Creates a NEW symbol, won't match dictionary key T_amb = self.parameters[sp.symbols("T_amb")] # KeyError! ``` **Why this fails:** Each call to `sp.symbols()` creates a new Python object with unique identity. The dictionary keys in `self.parameters` are the *original* symbol objects from `define_system()`, not new ones. **Solution 1: Store parameter values as attributes** (Recommended for beginners): ```python def define_system(self, T_amb_val=300.0, ...): # Store the numerical value self.T_amb_val = T_amb_val # Create symbol and use in parameters dict T_amb = sp.symbols('T_amb', positive=True) self.parameters = {T_amb: T_amb_val, ...} def setup_equilibria(self): # Use the stored value directly T_amb = self.T_amb_val # Works! ``` **Solution 2: Store symbol references** (Recommended for production): ```python def define_system(self, T_amb_val=300.0, ...): # Create symbol and store reference self.T_amb_sym = sp.symbols('T_amb', positive=True) self.parameters = {self.T_amb_sym: T_amb_val, ...} def setup_equilibria(self): # Use the stored symbol to look up value T_amb = self.parameters[self.T_amb_sym] # Works! ``` Our tutorial uses Solution 1 for simplicity and clarity. ::: ## Simulation and Visualization Now that we've defined our stochastic batch reactor, let's see it in action. This section demonstrates the complete workflow: inspecting the system, running simulations, and visualizing results. ::: {.callout-important} ## Integration Methods: `integrate()` vs `simulate()` ControlDESymulation provides two methods for time evolution: ### `integrate()` - Low-Level Interface **For stochastic systems, always use this method.** ```python result = system.integrate( x0=x0, u=controller, t_span=(0, 10), method='euler_maruyama', # Stochastic integrator dt=0.01, # Required for fixed-step methods seed=42, # Reproducibility n_paths=100 # Monte Carlo ensemble ) ``` **Stochastic-specific features:** - `seed`: Random number generator seed for reproducible noise - `n_paths`: Built-in Monte Carlo simulation (returns shape `(n_paths, T, nx)`) - Stochastic methods: `'euler_maruyama'`, `'milstein'`, `'heun'`, etc. - `dt`: Time step (required for most SDE methods) - `t_eval`: Optional array of specific output times **Fixed-Step vs Adaptive Methods:** Most SDE integrators are **fixed-step** (Euler-Maruyama, Milstein, Heun, etc.): ```python result = system.integrate( x0=x0, t_span=(0, 10), method='euler_maruyama', dt=0.01 # Required - sets uniform time grid ) # result['t'] has uniform spacing: [0.00, 0.01, 0.02, ..., 10.00] ``` Some advanced methods are **adaptive** (use `rtol`/`atol` instead of `dt`): ```python result = system.integrate( x0=x0, t_span=(0, 10), method='adaptive_sde', # Hypothetical adaptive method rtol=1e-6, atol=1e-8 ) # result['t'] may have non-uniform spacing ``` **Visualization tip:** For uniform time grids (most cases), plotting is straightforward. For adaptive grids, use `t_eval` to force specific output times: ```python result = system.integrate( x0=x0, t_span=(0, 10), method='adaptive_sde', t_eval=np.linspace(0, 10, 1000), # Force uniform grid rtol=1e-6 ) ``` **Returns:** `SDEIntegrationResult` with fields: - `t`: Time points `(T,)` - uniform for fixed-step, may vary for adaptive - `x`: States `(T, nx)` or `(n_paths, T, nx)` for Monte Carlo - `success`, `nfev`, `diffusion_evals`: Solver diagnostics ### `simulate()` - High-Level Interface **For deterministic systems only** (not recommended for SDEs). ```python result = system.simulate( x0=x0, controller=lambda x, t: -K @ x, t_span=(0, 10), dt=0.01 ) ``` **Key differences:** - Always returns regular time grid (uniform spacing `dt`) - Cleaner API for feedback controllers - No `seed` or `n_paths` support - Uses deterministic integrators (RK45, RK4, Euler) - Returns `SimulationResult` in time-major format `(T, nx)` **When to use each:** - **Stochastic systems (SDEs)** → `integrate()` with stochastic methods - **Deterministic systems (ODEs)** → Either method; `simulate()` is simpler - **Monte Carlo simulations** → `integrate()` with `n_paths > 1` For our stochastic batch reactor, we'll use `integrate()` exclusively. ::: ### Inspecting the System Before running simulations, let's examine what we've created: ```{python} # Display symbolic equations reactor = ContinuousStochasticBatchReactor( C_A0=1.0, # Initial concentration T0=350.0, # Initial temperature sigma_T=2.0 # Temperature noise intensity ) # Print the mathematical structure reactor.print_equations(simplify=True) ``` This shows: - State and control variables - System dimensions (nx=3 states, nu=1 control) - Symbolic drift dynamics with parameters substituted - Output map (if defined) We can also inspect system properties programmatically: ```{python} print(f"Number of states: {reactor.nx}") print(f"Number of controls: {reactor.nu}") print(f"State variables: {reactor.state_vars}") print(f"Parameters: {list(reactor.parameters.keys())}") ``` ::: {.callout-tip} ## Choosing an Integration Method For this tutorial, we use **Euler-Maruyama** (`method='euler_maruyama'`), the simplest fixed-step SDE method. It's a good starting point because: - **Fixed time step**: Requires `dt` parameter, produces uniform time grid - **Easy visualization**: Uniform spacing makes plotting straightforward - **Stable for additive noise**: Our reactor has additive noise (constant diffusion) - **Fast**: Minimal computation per step Other common methods: - **`milstein`**: Higher accuracy for multiplicative noise, still fixed-step - **`heun`**: Good balance of accuracy and stability, fixed-step - **Advanced methods**: May use adaptive stepping (see note below) ::: {.callout-note collapse="true"} ## Handling Adaptive SDE Methods Some advanced SDE solvers use adaptive time stepping for efficiency. If you use an adaptive method and need uniform output times (e.g., for certain visualizations or statistical analysis), use the `t_eval` parameter: ```python # Adaptive method with non-uniform internal steps result = system.integrate( x0=x0, t_span=(0, 10), method='adaptive_sde_method', # Hypothetical adaptive method rtol=1e-6, atol=1e-8, t_eval=np.linspace(0, 10, 1000) # Force uniform output grid ) ``` The solver will use adaptive steps internally for efficiency but interpolate to provide output at the exact times specified in `t_eval`. This gives you: - Efficiency of adaptive stepping during integration - Uniform time grid for downstream analysis and visualization - Control over output density **For this tutorial:** We use fixed-step Euler-Maruyama, so `dt` controls both internal steps and output times. No `t_eval` needed. ::: For production code with multiplicative noise, consider higher-order methods. For learning and visualization, Euler-Maruyama with `dt=0.01` to `dt=0.1` works well. ::: ### Equilibrium Points Our `setup_equilibria()` method created two equilibrium points. Let's see them: ```{python} # List all equilibria equilibrium_names = reactor.list_equilibria() print(f"Available equilibria: {equilibrium_names}") # Get the "complete conversion" equilibrium x_eq, u_eq = reactor.get_equilibrium('complete') print(f"\nComplete conversion equilibrium:") print(f" C_A = {x_eq[0]:.3f} mol/L") print(f" C_B = {x_eq[1]:.3f} mol/L") print(f" T = {x_eq[2]:.1f} K") print(f" Q = {u_eq[0]:.3f} W") # Verify it's actually an equilibrium (drift should be zero) drift = reactor(x_eq, u_eq, t=0) print(f"\nDrift at equilibrium: {drift}") print(f"Max drift magnitude: {np.abs(drift).max():.2e}") ``` **Key insight:** For stochastic systems, equilibria are fixed points of the **drift** only. The diffusion causes trajectories to fluctuate around these points with variance growing as $\sqrt{t}$. ## Summary: Building Your System To define a symbolic system in ControlDESymulation: **Minimum viable system:** ```python def define_system(self): # 1. Define symbolic variables x = sp.symbols('x') u = sp.symbols('u') # 2. Set required fields self.state_vars = [x] self.control_vars = [u] # or [] for autonomous self.parameters = {} # or {param: value} self._f_sym = sp.Matrix([...]) # Your dynamics # 3. For stochastic systems, add: # self.diffusion_expr = sp.Matrix([...]) ``` **Autonomous system (no control inputs):** ```python def define_system(self, damping=0.1, stiffness=1.0): # 1. Define symbolic variables (no control symbols needed) x, v = sp.symbols('x v') b, k = sp.symbols('b k', positive=True) # 2. Set required fields self.state_vars = [x, v] self.control_vars = [] # Empty list for autonomous self.parameters = {b: damping, k: stiffness} # 3. Dynamics depending only on states self._f_sym = sp.Matrix([ v, # dx/dt = v -k*x - b*v # dv/dt = -kx - bv (damped oscillator) ]) ``` **Pure diffusion system (stochastic, zero drift):** ```python def define_system(self, sigma_x=1.0, sigma_y=1.0, sigma_z=1.0): # 1. Define symbolic variables x, y, z = sp.symbols('x y z', real=True) sx, sy, sz = sp.symbols('sigma_x sigma_y sigma_z', positive=True) # 2. Set required fields self.state_vars = [x, y, z] self.control_vars = [] # Typically autonomous self.parameters = {sx: sigma_x, sy: sigma_y, sz: sigma_z} # 3. ZERO DRIFT: Pure diffusion has no deterministic dynamics self._f_sym = sp.zeros(3, 1) # Must be zero vector # 4. DIFFUSION: Stochastic dynamics only self.diffusion_expr = sp.Matrix([ [sx, 0, 0], # Independent noise in x [0, sy, 0], # Independent noise in y [0, 0, sz] # Independent noise in z ]) self.sde_type = "ito" ``` **Well-documented system (recommended):** ```python def define_system(self, param1=1.0, param2=2.0): # 1. Define all symbolic variables (kept as LOCAL variables) x, y = sp.symbols('x y') # State symbols - local only u = sp.symbols('u') # Control symbol - local only p1, p2 = sp.symbols('p1 p2', positive=True) # Parameter symbols - local only # 2. Set required fields self.state_vars = [x, y] # Reference to local symbols self.control_vars = [u] # Reference to local symbols self.parameters = {p1: param1, p2: param2} # Symbol→value mapping self._f_sym = sp.Matrix([...]) # 3. Explicitly set optional fields for clarity self.output_vars = [] # No computed outputs self._h_sym = sp.Matrix([x, y]) # Observe full state self.order = 1 # First-order form # 4. For stochastic systems: # self.diffusion_expr = sp.Matrix([...]) # self.sde_type = "ito" # or "stratonovich" # 5. Store parameter VALUES as attributes (if needed in setup_equilibria) # DON'T store symbols unless parameters need to change post-initialization # symbols are already in self.state_vars, self.parameters # self.param1_val = param1 # Store VALUE for later use # self.param2_val = param2 # Store VALUE for later use def setup_equilibria(self): # Optional: Add parameter-dependent equilibria # Access stored values: x_eq_value = self.param1_val self.add_equilibrium("my_equilibrium", x_eq=..., u_eq=...) ``` **Which template should you use?** - **Minimum viable**: Quick prototyping, testing, simple demos - **Autonomous system**: No control inputs (free oscillators, population models, unforced reactions) - **Pure diffusion**: Stochastic processes with zero drift (Brownian motion, pure noise) - **Well-documented**: Production code, complex systems, collaborative projects (recommended) The reactor example in this tutorial follows the well-documented approach with both drift and diffusion. ## What You've Learned This tutorial covered the **anatomy of a symbolic system** through a stochastic batch reactor example: ### System Definition - How to subclass `ContinuousStochasticSystem` - Required fields: `state_vars`, `control_vars`, `parameters`, `_f_sym`, `diffusion_expr` - Optional fields: `output_vars`, `_h_sym`, `order`, `sde_type` - Special cases: Autonomous systems and pure diffusion systems - Parameter-dependent equilibria via `setup_equilibria()` - **What to store as attributes vs local variables**: - Symbols (state, control, parameters) → Keep local - Parameter values (used in `setup_equilibria`) → Store as `self.param_val` - Initial conditions and custom data → Store as attributes ### System Introspection - `print_equations()`: Display symbolic dynamics - `list_equilibria()`: See all equilibrium points - `get_equilibrium()`: Access equilibrium states and controls - Programmatic access: `nx`, `nu`, `state_vars`, `parameters` ### Critical Distinctions - **Drift vs Diffusion**: Deterministic dynamics vs stochastic forcing - **Itô vs Stratonovich**: Different stochastic calculus interpretations - **Equilibria in SDEs**: Fixed points of drift only (stochastic fluctuations persist) - **integrate() vs simulate()**: Use `integrate()` for stochastic systems - **Fixed-step vs adaptive**: Most SDE methods use fixed `dt` (uniform time grid), some advanced methods adapt step size (use `t_eval` for uniform output) - **Parameter access**: Store parameter values as attributes (`self.param_val`) for use in `setup_equilibria()`, don't try to look up via `self.parameters[sp.symbols("param")]` (creates new symbol, causes KeyError). Store symbol references as attributes only if parameter values need to change post-initialization. ## Next Steps **[Tutorial 2: Simulation and Visualization](batch_reactor_simulation.qmd)** - Single trajectory simulation with `integrate()` - Monte Carlo ensembles with `n_paths` - Time series visualization with `reactor.plotter` - Phase portraits and state-space analysis - Ensemble statistics and confidence bands

System Overview

Noise Structure

What This Tutorial Covers

System Definition

SymPy Symbols: Usually Local, Sometimes Stored

Parameter Values: Store When Needed Later

Custom Fields: Always Store

Decision Tree

Required Specification

Dynamics for Autonomous Systems

Calling Conventions

Converting Between Controlled and Autonomous

Time-Varying vs Autonomous

Practical Examples

Mathematical Form

Required Specification

Example: 3D Brownian Motion

When to Use Pure Diffusion

Physical Interpretation

Integration Considerations

Drift + Diffusion vs Pure Diffusion

Practical Applications

Required vs Optional System Definition Elements

self.output_vars (default: [])

self._h_sym (default: full state vector)

self.order (default: 1)

self.sde_type (default: "ito")

setup_equilibria()

Custom Methods and Fields

Understanding System Order in Detail

First-Order State-Space Form (order=1, default)

Higher-Order Form (order=n)

When to Use Each Form

Requirements

This Tutorial

Equilibrium Setup Example

Simulation and Visualization

integrate() - Low-Level Interface

simulate() - High-Level Interface

Inspecting the System

Equilibrium Points

Summary: Building Your System

What You’ve Learned

System Definition

System Introspection

Critical Distinctions

Next Steps

`self.output_vars` (default: `[]`)

`self._h_sym` (default: full state vector)

`self.order` (default: `1`)

`self.sde_type` (default: `"ito"`)

`setup_equilibria()`

`integrate()` - Low-Level Interface

`simulate()` - High-Level Interface