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¶
[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 pathlib import Path
from pathsim import Simulation, Connection
from pathsim.blocks import ModelExchangeFMU, Scope
from pathsim.solvers import RKDP54, ESDIRK43
FMU Path¶
[2]:
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.
[3]:
# 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:
[4]:
# Display FMU metadata (accessed via fmu_wrapper)
md = fmu.fmu_wrapper.model_description
print(f"Model Name: {md.modelName}")
print(f"FMI Version: {fmu.fmu_wrapper.fmi_version}")
print(f"Description: {md.description}")
print(f"Generation Tool: {md.generationTool}")
print(f"\nContinuous states: {fmu.fmu_wrapper.n_states}")
print(f"Event indicators: {fmu.fmu_wrapper.n_event_indicators}")
print(f"Outputs: {len(fmu.fmu_wrapper.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.
[5]:
# 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)
07:41:21 - INFO - LOGGING (log: True)
07:41:21 - INFO - BLOCKS (total: 2, dynamic: 1, static: 1, eventful: 0)
07:41:21 - INFO - GRAPH (nodes: 2, edges: 2, alg. depth: 1, loop depth: 0, runtime: 0.130ms)
07:41:21 - INFO - STARTING -> TRANSIENT (Duration: 30.00s)
07:41:21 - INFO - -------------------- 1% | 0.0s<0.2s | 1873.2 it/s
07:41:21 - INFO - ####---------------- 20% | 0.0s<0.1s | 2279.5 it/s
07:41:21 - INFO - ########------------ 40% | 0.1s<0.1s | 2300.4 it/s
07:41:21 - INFO - ############-------- 60% | 0.1s<0.1s | 2202.1 it/s
07:41:21 - INFO - ################---- 80% | 0.2s<0.0s | 2321.9 it/s
07:41:21 - INFO - #################### 100% | 0.2s<--:-- | 2237.9 it/s
07:41:21 - INFO - FINISHED -> TRANSIENT (total steps: 415, successful: 351, runtime: 195.01 ms)
[5]:
{'total_steps': 415, 'successful_steps': 351, 'runtime_ms': 195.00624600004812}
Results¶
[6]:
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.
[7]:
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)
07:41:21 - INFO - LOGGING (log: True)
07:41:21 - INFO - BLOCKS (total: 2, dynamic: 1, static: 1, eventful: 0)
07:41:21 - INFO - GRAPH (nodes: 2, edges: 2, alg. depth: 1, loop depth: 0, runtime: 0.138ms)
07:41:21 - INFO - STARTING -> TRANSIENT (Duration: 1500.00s)
07:41:21 - INFO - FINISHED -> TRANSIENT (total steps: 0, successful: 0, runtime: 5.01 ms)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
ValueError: On entry to DLASCL parameter number 4 had an illegal value
The above exception was the direct cause of the following exception:
SystemError Traceback (most recent call last)
Cell In[7], line 22
18 tolerance_fpi=1e-8,
19 log=True
20 )
21
---> 22 sim_stiff.run(1500.0)
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.21.0/lib/python3.13/site-packages/pathsim/simulation.py:1723, in Simulation.run(self, duration, reset, adaptive)
1721 #enter tracker context and consume the run loop
1722 with tracker:
-> 1723 for _ in self._run_loop(duration, reset, adaptive, tracker=tracker):
1724 pass
1726 return tracker.stats
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.21.0/lib/python3.13/site-packages/pathsim/simulation.py:1644, in Simulation._run_loop(self, duration, reset, adaptive, tracker)
1640 #main simulation loop
1641 while self.time < end_time and self._active:
1642
1643 #advance the simulation by one (effective) timestep '_dt'
-> 1644 success, error_norm, scale, *_ = self.timestep(
1645 dt=_dt,
1646 adaptive=_adaptive
1647 )
1649 #perform adaptive rescale
1650 if _adaptive:
1651
1652 #if no error estimate and rescale -> back to default timestep
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.21.0/lib/python3.13/site-packages/pathsim/simulation.py:1469, in Simulation.timestep(self, dt, adaptive)
1465 for time_stage in self.engine.stages(self.time, dt):
1467 if is_implicit:
1468 #implicit: solve update equation (contains _update internally)
-> 1469 success, evals, solver_its = self._solve(time_stage, dt)
1470 total_evals += evals
1471 total_solver_its += solver_its
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.21.0/lib/python3.13/site-packages/pathsim/simulation.py:1113, in Simulation._solve(self, t, dt)
1111 if not block:
1112 continue
-> 1113 self._solve_tracker.record(block, block.solve(t, dt))
1115 #check for convergence
1116 if self._solve_tracker.converged(self.tolerance_fpi):
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.21.0/lib/python3.13/site-packages/pathsim/blocks/dynsys.py:131, in DynamicalSystem.solve(self, t, dt)
129 x, u = self.engine.state, self.inputs.to_array()
130 f, J = self.op_dyn(x, u, t), self.op_dyn.jac_x(x, u, t)
--> 131 return self.engine.solve(f, J, dt)
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.21.0/lib/python3.13/site-packages/pathsim/solvers/_rungekutta.py:312, in DiagonallyImplicitRungeKutta.solve(self, f, J, dt)
309 b = self.BT[self.stage][self.stage]
311 #optimizer step with block local jacobian
--> 312 self.x, err = self.opt.step(self.x, x_0 + dt * slope, dt * b * J)
314 else:
315 #optimizer step (pure)
316 self.x, err = self.opt.step(self.x, x_0 + dt * slope, None)
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.21.0/lib/python3.13/site-packages/pathsim/optim/anderson.py:365, in NewtonAnderson.step(self, x, g, jac)
362 _x, res_norm = self._newton(x, g, jac)
364 #anderson step with no residual
--> 365 y, _ = super().step(_x, g)
367 return y, res_norm
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.21.0/lib/python3.13/site-packages/pathsim/optim/anderson.py:213, in Anderson.step(self, x, g)
210 return _x - _res[0] * np.dot(dR_flat, dX_flat) / dR2, abs(_res[0])
212 #compute coefficients from least squares problem
--> 213 C, *_ = np.linalg.lstsq(dR.T, _res, rcond=None)
215 #new solution and residual norm
216 return _x - C @ dX, np.linalg.norm(_res)
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.21.0/lib/python3.13/site-packages/numpy/linalg/linalg.py:2326, in lstsq(a, b, rcond)
2323 if n_rhs == 0:
2324 # lapack can't handle n_rhs = 0 - so allocate the array one larger in that axis
2325 b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype)
-> 2326 x, resids, rank, s = gufunc(a, b, rcond, signature=signature, extobj=extobj)
2327 if m == 0:
2328 x[...] = 0
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.21.0/lib/python3.13/site-packages/numpy/linalg/linalg.py:124, in _raise_linalgerror_lstsq(err, flag)
123 def _raise_linalgerror_lstsq(err, flag):
--> 124 raise LinAlgError("SVD did not converge in Linear Least Squares")
SystemError: <class 'numpy.linalg.LinAlgError'> returned a result with an exception set
[8]:
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