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)
14:02:24 - INFO - LOGGING (log: True)
14:02:24 - INFO - BLOCKS (total: 2, dynamic: 1, static: 1, eventful: 0)
14:02:24 - INFO - GRAPH (nodes: 2, edges: 2, alg. depth: 1, loop depth: 0, runtime: 0.125ms)
14:02:24 - INFO - STARTING -> TRANSIENT (Duration: 30.00s)
14:02:24 - INFO - -------------------- 1% | 0.0s<0.1s | 3021.5 it/s
14:02:24 - INFO - ####---------------- 20% | 0.0s<0.1s | 3734.1 it/s
14:02:24 - INFO - ########------------ 40% | 0.0s<0.0s | 3761.0 it/s
14:02:24 - INFO - ############-------- 60% | 0.1s<0.0s | 3775.6 it/s
14:02:24 - INFO - ################---- 80% | 0.1s<0.0s | 3698.8 it/s
14:02:24 - INFO - #################### 100% | 0.1s<--:-- | 3811.3 it/s
14:02:24 - INFO - FINISHED -> TRANSIENT (total steps: 415, successful: 351, runtime: 115.59 ms)
[5]:
{'total_steps': 415, 'successful_steps': 351, 'runtime_ms': 115.58894199970382}
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)
14:02:24 - INFO - LOGGING (log: True)
14:02:24 - INFO - BLOCKS (total: 2, dynamic: 1, static: 1, eventful: 0)
14:02:24 - INFO - GRAPH (nodes: 2, edges: 2, alg. depth: 1, loop depth: 0, runtime: 0.089ms)
14:02:24 - INFO - STARTING -> TRANSIENT (Duration: 1500.00s)
14:02:24 - INFO - FINISHED -> TRANSIENT (total steps: 0, successful: 0, runtime: 2.61 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
9 sco_stiff = Scope(labels=[r"$x_0$", r"$x_1$"])
11 sim_stiff = Simulation(
12 [fmu_stiff, sco_stiff],
13 [Connection(fmu_stiff[0], sco_stiff[0]), Connection(fmu_stiff[1], sco_stiff[1])],
(...) 19 log=True
20 )
---> 22 sim_stiff.run(1500.0)
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.15.2/lib/python3.13/site-packages/pathsim/simulation.py:1631, in Simulation.run(self, duration, reset, adaptive)
1629 #enter tracker context and consume the run loop
1630 with tracker:
-> 1631 for _ in self._run_loop(duration, reset, adaptive, tracker=tracker):
1632 pass
1634 return tracker.stats
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.15.2/lib/python3.13/site-packages/pathsim/simulation.py:1552, in Simulation._run_loop(self, duration, reset, adaptive, tracker)
1548 #main simulation loop
1549 while self.time < end_time and self._active:
1550
1551 #advance the simulation by one (effective) timestep '_dt'
-> 1552 success, error_norm, scale, *_ = self.timestep(
1553 dt=_dt,
1554 adaptive=_adaptive
1555 )
1557 #perform adaptive rescale
1558 if _adaptive:
1559
1560 #if no error estimate and rescale -> back to default timestep
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.15.2/lib/python3.13/site-packages/pathsim/simulation.py:1447, in Simulation.timestep(self, dt, adaptive)
1445 return self.timestep_adaptive_explicit(dt)
1446 else:
-> 1447 return self.timestep_adaptive_implicit(dt)
1448 else:
1449 if self.engine.is_explicit:
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.15.2/lib/python3.13/site-packages/pathsim/simulation.py:1364, in Simulation.timestep_adaptive_implicit(self, dt)
1358 if self._blocks_dyn:
1359
1360 #iterate explicit solver stages with evaluation time (generator)
1361 for time_stage in self.engine.stages(self.time, dt):
1362
1363 #solve implicit update equation and get iteration count
-> 1364 success, evals, solver_its = self._solve(time_stage, dt)
1366 #count solver iterations and function evaluations
1367 total_solver_its += solver_its
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.15.2/lib/python3.13/site-packages/pathsim/simulation.py:834, in Simulation._solve(self, t, dt)
831 continue
833 #advance solution (internal optimizer)
--> 834 error = block.solve(t, dt)
835 if error > max_error:
836 max_error = error
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.15.2/lib/python3.13/site-packages/pathsim/blocks/dynsys.py:151, in DynamicalSystem.solve(self, t, dt)
149 x, u = self.engine.get(), self.inputs.to_array()
150 f, J = self.op_dyn(x, u, t), self.op_dyn.jac_x(x, u, t)
--> 151 return self.engine.solve(f, J, dt)
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.15.2/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.15.2/lib/python3.13/site-packages/pathsim/optim/anderson.py:317, in NewtonAnderson.step(self, x, g, jac)
314 _x, res_norm = self._newton(x, g, jac)
316 #anderson step with no residual
--> 317 y, _ = super().step(_x, g)
319 return y, res_norm
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.15.2/lib/python3.13/site-packages/pathsim/optim/anderson.py:187, in Anderson.step(self, x, g)
184 return _x - _res * np.dot(dR, dX) / dR2, abs(_res)
186 #compute coefficients from least squares problem
--> 187 C, *_ = np.linalg.lstsq(dR.T, _res, rcond=None)
189 #new solution and residual norm
190 return _x - C @ dX, np.linalg.norm(_res)
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.15.2/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.15.2/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