Nested Subsystems

Demonstration of hierarchical modeling using nested subsystems for a Van der Pol oscillator.

You can also find this example as a single file in the GitHub repository.

Why Use Subsystems?

Subsystems allow you to:

  • Organize complex systems into logical modules

  • Reuse components across different models

  • Abstract implementation details

  • Scale to large systems with many components

  • Debug and test individual modules separately

The Van der Pol Oscillator Revisited

The stiff Van der Pol oscillator is described by:

\[\frac{dx_1}{dt} = x_2\]
\[\frac{dx_2}{dt} = \mu(1 - x_1^2)x_2 - x_1\]

With \(\mu = 1000\) (very stiff!)

Hierarchical Structure

This example demonstrates hierarchical modeling using Subsystem and Interface blocks for modular system design.

[1]:

import numpy as np
import matplotlib.pyplot as plt

# Apply PathSim docs matplotlib style for consistent, theme-friendly figures
plt.style.use('../pathsim_docs.mplstyle')

from pathsim import Simulation, Connection, Interface, Subsystem
from pathsim.blocks import Integrator, Scope, Function, Multiplier, Adder, Amplifier, Constant, Pow
from pathsim.solvers import ESDIRK43

System Parameters

[2]:

# Initial conditions
x1_0 = 2.0
x2_0 = 0.0

# Van der Pol parameter (high stiffness!)
mu = 1000

# Simulation timestep
dt = 0.01

Level 1: ODE Function Subsystem

First, we create a subsystem that computes \(f(x_1, x_2) = \mu(1 - x_1^2)x_2 - x_1\)

This subsystem:

  • Takes two inputs: \(x_1\) and \(x_2\)

  • Returns one output: the computed derivative

  • Is self-contained and reusable

vanderpol function as a subsystem

The Interface block defines the subsystem’s inputs and outputs and this is how it looks like in pathsim:

[3]:

# Interface for the ODE function subsystem
# Input 0: x1
# Input 1: x2
# Output: μ(1 - x1²)x2 - x1

In = Interface()
M1 = Multiplier()    # For x2 * (1 - x1²)
C1 = Constant(1)     # The constant 1
A1 = Amplifier(mu)   # Multiply by μ
P1 = Adder("+-")     # Sum: μ(1 - x1²)x2 - x1
P2 = Adder("+-")     # Compute: 1 - x1²
S1 = Pow(2)          # Square x1

fn_blocks = [In, M1, C1, A1, P1, P2, S1]

fn_connections = [
    Connection(In[0], S1),         # x1 → x1²
    Connection(In[1], M1[0]),      # x2 → multiplier
    Connection(S1, P2[1]),         # x1² to adder (-)
    Connection(C1, P2[0]),         # 1 to adder
    Connection(P2, M1[1]),         # (1 - x1²) to multiplier
    Connection(M1, A1),            # x2(1 - x1²) → multiply by μ
    Connection(A1, P1[0]),         # μ(...) to final adder
    Connection(In[0], P1[1]),      # x1 to final adder (-)
    Connection(P1, In)             # Result to output
]

# Create the ODE function subsystem
Fn = Subsystem(fn_blocks, fn_connections)

Level 2: Van der Pol Subsystem

Now we create a subsystem that contains:

  • Two integrators (for \(x_1\) and \(x_2\))

  • The ODE function subsystem we just created

vanderpol ODE as a subsystem

This implements the complete Van der Pol ODE system:

[4]:

# Interface for VDP subsystem (no inputs, two outputs)
If = Interface()
I1 = Integrator(x1_0)  # Integrator for x1
I2 = Integrator(x2_0)  # Integrator for x2

vdp_blocks = [If, I1, I2, Fn]

vdp_connections = [
    Connection(I2, I1, Fn[1], If[1]),  # x2 → I1, Fn, and output
    Connection(I1, Fn, If),            # x1 → Fn and output[0]
    Connection(Fn, I2)                 # dx2/dt → I2
]

# Create the Van der Pol subsystem
VDP = Subsystem(vdp_blocks, vdp_connections)

Level 3: Top-Level System

Finally, we create the top-level system that contains:

  • The VDP subsystem

  • A Scope for visualization

top level system view

At this level, the VDP subsystem looks like a simple block with two outputs, hiding all its internal complexity

[5]:

# Top-level system
Sco = Scope(labels=[r"$x_1(t)$", r"$x_2(t)$"])

blocks = [VDP, Sco]

connections = [
    Connection(VDP, Sco),       # x1 to scope
    # Connection(VDP[1], Sco[1])  # x2 to scope
]

Simulation Setup

We use a stiff solver (ESDIRK43) because :math:`mu = 1000makes this a very stiff system.

[6]:

# Initialize simulation
Sim = Simulation(
    blocks,
    connections,
    dt=dt,
    log=True,
    Solver=ESDIRK43,
    tolerance_lte_abs=1e-6,
    tolerance_lte_rel=1e-4,
    tolerance_fpi=1e-7
)

# Run simulation (note: 2*mu for long-term dynamics)
Sim.run(2*mu)

07:09:33 - INFO - LOGGING (log: True)
07:09:33 - INFO - BLOCKS (total: 2, dynamic: 1, static: 1, eventful: 0)
07:09:33 - INFO - GRAPH (nodes: 2, edges: 1, alg. depth: 1, loop depth: 0, runtime: 0.022ms)
07:09:33 - INFO - STARTING -> TRANSIENT (Duration: 2000.00s)
07:09:33 - INFO - --------------------   1% | 0.4s<11.4s | 45.5 it/s
07:09:33 - WARNING - implicit solver not converged, reverting step at t=557.471131
  Subsystem: 1.19e-01
07:09:33 - INFO - #####---------------  27% | 0.7s<1.1s | 27.6 it/s
07:09:34 - WARNING - implicit solver not converged, reverting step at t=676.479189
  Subsystem: 8.00e-03
07:09:34 - WARNING - implicit solver not converged, reverting step at t=744.641850
  Subsystem: 5.03e-02
07:09:34 - WARNING - implicit solver not converged, reverting step at t=716.038014
  Subsystem: 1.11e-04
07:09:34 - WARNING - implicit solver not converged, reverting step at t=754.836313
  Subsystem: 3.80e-04
07:09:34 - WARNING - implicit solver not converged, reverting step at t=784.117711
  Subsystem: 3.89e-02
07:09:34 - WARNING - implicit solver not converged, reverting step at t=797.308454
  Subsystem: 1.49e-02
07:09:34 - WARNING - implicit solver not converged, reverting step at t=800.647586
  Subsystem: 2.40e-05
07:09:34 - INFO - ########------------  40% | 1.2s<8.0s | 29.4 it/s
07:09:34 - WARNING - implicit solver not converged, reverting step at t=803.643644
  Subsystem: 8.47e-07
07:09:34 - WARNING - implicit solver not converged, reverting step at t=805.302856
  Subsystem: 3.50e-06
07:09:34 - WARNING - implicit solver not converged, reverting step at t=806.501125
  Subsystem: 6.08e-04
07:09:36 - INFO - ########------------  41% | 2.9s<6.8s | 67.4 it/s
07:09:36 - WARNING - implicit solver not converged, reverting step at t=1250.693446
  Subsystem: 1.79e-03
07:09:36 - INFO - ############--------  62% | 3.1s<0.6s | 34.1 it/s
07:09:36 - WARNING - implicit solver not converged, reverting step at t=1412.300094
  Subsystem: 2.35e-03
07:09:36 - WARNING - implicit solver not converged, reverting step at t=1331.496770
  Subsystem: 3.14e-06
07:09:36 - WARNING - implicit solver not converged, reverting step at t=1443.543101
  Subsystem: 3.00e-05
07:09:36 - WARNING - implicit solver not converged, reverting step at t=1525.065843
  Subsystem: 7.51e-03
07:09:36 - WARNING - implicit solver not converged, reverting step at t=1484.304472
  Subsystem: 3.73e-05
07:09:36 - WARNING - implicit solver not converged, reverting step at t=1541.546358
  Subsystem: 9.09e-03
07:09:36 - WARNING - implicit solver not converged, reverting step at t=1580.432813
  Subsystem: 2.23e-03
07:09:36 - WARNING - implicit solver not converged, reverting step at t=1580.432813
  Subsystem: 3.48e-04
07:09:36 - WARNING - implicit solver not converged, reverting step at t=1585.468308
  Subsystem: 7.52e-03
07:09:37 - WARNING - implicit solver not converged, reverting step at t=1600.952536
  Subsystem: 5.22e-02
07:09:37 - WARNING - implicit solver not converged, reverting step at t=1594.454722
  Subsystem: 1.03e-05
07:09:37 - WARNING - implicit solver not converged, reverting step at t=1611.057481
  Subsystem: 2.58e+00
07:09:37 - INFO - ################----  80% | 3.9s<3.4s | 21.6 it/s
07:09:37 - WARNING - implicit solver not converged, reverting step at t=1607.880672
  Subsystem: 3.77e-05
07:09:37 - WARNING - implicit solver not converged, reverting step at t=1607.880672
  Subsystem: 4.21e-04
07:09:37 - WARNING - implicit solver not converged, reverting step at t=1610.569142
  Subsystem: 4.27e-06
07:09:37 - WARNING - implicit solver not converged, reverting step at t=1612.041876
  Subsystem: 1.60e-06
07:09:37 - WARNING - implicit solver not converged, reverting step at t=1613.632269
  Subsystem: 1.52e-03
07:09:39 - INFO - ################----  81% | 5.7s<3.8s | 79.4 it/s
07:09:39 - INFO - #################### 100% | 5.9s<--:-- | 38.1 it/s
07:09:39 - INFO - FINISHED -> TRANSIENT (total steps: 363, successful: 246, runtime: 5894.86 ms)

[6]:

{'total_steps': 363, 'successful_steps': 246, 'runtime_ms': 5894.859768999595}

Results: Time Series

The Van der Pol oscillator with \(\mu = 1000\) exhibits relaxation oscillations - fast transitions between slow phases. This requires a stiff solver to handle efficiently.

[7]:

Sco.plot(".-", lw=1.5)
plt.show()
../_images/examples_nested_subsystems_19_0.svg

Subsystem Benefits Demonstrated

This example shows several advantages of subsystems:

  1. Modularity: The ODE function is completely separate from the integration

  2. Reusability: The Fn subsystem could be used in other models

  3. Clarity: The top level is clean - just VDP and Scope

  4. Debugging: Each subsystem can be tested independently

  5. Abstraction: Inner complexity is hidden from higher levels

Comparison with ODE Block

Compare this hierarchical approach with using a single ODE block:

# Alternative: Using ODE block (simpler but less modular)
def vdp_ode(x, u, t):
    return np.array([x[1], mu*(1 - x[0]**2)*x[1] - x[0]])

VDP = ODE(vdp_ode, np.array([x1_0, x2_0]))

Both approaches work! Use subsystems when:

  • You need modularity and reusability

  • The system is complex with many components

  • You want to visualize internal signals

  • You’re building block diagram models