FMU ME: Van der Pol

This example demonstrates Model Exchange FMU integration with a nonlinear oscillator. The Van der Pol equation exhibits self-sustained oscillations:

\[\ddot{x} - \mu(1 - x^2)\dot{x} + x = 0\]

where \(\mu > 0\) controls nonlinearity. As a first-order system:

\[\frac{dx_0}{dt} = x_1\]

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

You can also find the FMU integration tests in the GitHub repository.

This example demonstrates the ModelExchangeFMU block with purely continuous-time dynamics (no events), showcasing PathSim’s adaptive integration for nonlinear systems.

Import and Setup

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 pathlib import Path
from pathsim import Simulation, Connection
from pathsim.blocks import ModelExchangeFMU, Scope
from pathsim.solvers import RKDP54, ESDIRK43

FMU Path

notebook_dir = Path().resolve()
fmu_path = notebook_dir / "data" / "VanDerPol_ME.fmu"

# Verify FMU exists
if not fmu_path.exists():
    raise FileNotFoundError(f"FMU file not found at {fmu_path}")

System Definition

We simulate with \(\mu = 1.0\) for clear nonlinear oscillations:

The ModelExchangeFMU exposes states $(x_0, x_1)$ and provides derivatives to PathSim’s solvers.

# Create Model Exchange FMU
fmu = ModelExchangeFMU(
    fmu_path=str(fmu_path),
    instance_name="vanderpol",
    start_values={
        "mu": 1.0,     # nonlinearity parameter
        "x0": 2.0,     # initial position
        "x1": 0.0,     # initial velocity
    },
    tolerance=1e-8,
    verbose=False
)

# Scope to record states
sco = Scope(labels=[r"$x_0$", r"$x_1$"])

blocks = [fmu, sco]

# Connections
connections = [
    Connection(fmu[0], sco[0]),
    Connection(fmu[1], sco[1]),
]

Display FMU metadata:

print(f"Model Name: {fmu.model_name}")
print(f"FMI Version: {fmu.fmi_version}")
print(f"Description: {fmu.description}")
print(f"Generation Tool: {fmu.generation_tool}")
print(f"\nContinuous states: {fmu.n_states}")
print(f"Event indicators: {fmu.n_event_indicators}")
print(f"Outputs: {len(fmu._output_refs)}")

Model Name: Van der Pol oscillator
FMI Version: 2.0
Description: This model implements the van der Pol oscillator
Generation Tool: Reference FMUs (v0.0.39)

Continuous states: 2
Event indicators: 0
Outputs: 2

Simulation Setup

We use a high-order adaptive solver (RKDP54) for accurate nonlinear integration.

# Initialize simulation
sim = Simulation(
    blocks,
    connections,
    dt=0.1,
    dt_max=0.1,
    Solver=RKDP54,
    tolerance_lte_rel=1e-6,
    tolerance_lte_abs=1e-9,
    log=True
)

# Run simulation
sim.run(30.0)

22:36:39 - INFO - LOGGING (log: True)
22:36:39 - INFO - BLOCKS (total: 2, dynamic: 1, static: 1, eventful: 0)
22:36:39 - INFO - GRAPH (nodes: 2, edges: 2, alg. depth: 1, loop depth: 0, runtime: 0.220ms)
22:36:39 - INFO - STARTING -> TRANSIENT (Duration: 30.00s)
22:36:39 - INFO - --------------------   1% | 0.0s<0.2s | 2976.8 it/s
22:36:39 - INFO - ####----------------  20% | 0.0s<0.1s | 3794.0 it/s
22:36:39 - INFO - ########------------  40% | 0.0s<0.0s | 3744.9 it/s
22:36:39 - INFO - ############--------  60% | 0.1s<0.0s | 3697.6 it/s
22:36:39 - INFO - ################----  80% | 0.1s<0.0s | 3722.5 it/s
22:36:39 - INFO - #################### 100% | 0.1s<--:-- | 3736.6 it/s
22:36:39 - INFO - FINISHED -> TRANSIENT (total steps: 415, successful: 351, runtime: 115.33 ms)

{'total_steps': 415, 'successful_steps': 351, 'runtime_ms': 115.32791200079373}

Results

sco.plot()
plt.show()

Relaxation Oscillations

For very large \(\mu\), the system exhibits relaxation oscillations with extreme stiffness. We use \(\mu = 500\) to demonstrate PathSim’s capability with severe stiffness:

Stiff systems require implicit solvers with large stability regions. We use ESDIRK43, a 4th-order implicit solver designed for stiff problems.

fmu_stiff = ModelExchangeFMU(
    fmu_path=str(fmu_path),
    instance_name="vdp_stiff",
    start_values={"mu": 500.0, "x0": 2.0, "x1": 0.0},
    tolerance=1e-8,
    verbose=False
)

sco_stiff = Scope(labels=[r"$x_0$", r"$x_1$"])

sim_stiff = Simulation(
    [fmu_stiff, sco_stiff],
    [Connection(fmu_stiff[0], sco_stiff[0]), Connection(fmu_stiff[1], sco_stiff[1])],
    dt=0.1,
    Solver=ESDIRK43,          # implicit solver for stiff systems
    tolerance_lte_rel=1e-4,
    tolerance_lte_abs=1e-6,
    tolerance_fpi=1e-8,
    log=True
)

sim_stiff.run(1500.0)

22:36:39 - INFO - LOGGING (log: True)
22:36:39 - INFO - BLOCKS (total: 2, dynamic: 1, static: 1, eventful: 0)
22:36:39 - INFO - GRAPH (nodes: 2, edges: 2, alg. depth: 1, loop depth: 0, runtime: 0.161ms)
22:36:39 - INFO - STARTING -> TRANSIENT (Duration: 1500.00s)
22:36:39 - INFO - #-------------------   9% | 0.1s<0.2s | 224.3 it/s
22:36:39 - INFO - ####----------------  21% | 0.2s<0.2s | 134.5 it/s
22:36:40 - INFO - #########-----------  45% | 1.0s<0.1s | 144.8 it/s
22:36:41 - INFO - #############-------  66% | 1.7s<0.1s | 190.1 it/s
22:36:41 - INFO - ################----  80% | 1.9s<0.2s | 110.2 it/s
22:36:42 - INFO - #################### 100% | 2.6s<--:-- | 179.5 it/s
22:36:42 - INFO - FINISHED -> TRANSIENT (total steps: 372, successful: 237, runtime: 2565.11 ms)

{'total_steps': 372, 'successful_steps': 237, 'runtime_ms': 2565.1059330011776}

time_stiff, (x0_stiff, x1_stiff) = sco_stiff.read()

fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True, dpi=130)

axes[0].plot(time_stiff, x0_stiff, lw=2)
axes[0].set_ylabel(r'$x_0$')
axes[0].set_title(r'Relaxation Oscillations ($\mu = 500$)')
axes[0].grid(True, alpha=0.3)

axes[1].plot(time_stiff, x1_stiff, lw=2, color='orange')
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel(r'$x_1$')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Key Features

• Seamless integration of nonlinear dynamics

• Adaptive timestepping for varying dynamics

• Parameter sweep capabilities

• Handles moderately stiff systems