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]),
]
/tmp/ipykernel_5265/203168073.py:2: DeprecationWarning: 'ModelExchangeFMU' is deprecated and will be removed in version 1.0.0. Use 'pathsim_fmi.ModelExchangeFMU' instead. This block has moved to the pathsim-fmi package: pip install pathsim-fmi
fmu = ModelExchangeFMU(
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)
11:40:08 - INFO - LOGGING (log: True)
11:40:08 - INFO - BLOCKS (total: 2, dynamic: 1, static: 1, eventful: 0)
11:40:08 - INFO - GRAPH (nodes: 2, edges: 2, alg. depth: 1, loop depth: 0, runtime: 0.128ms)
11:40:08 - INFO - STARTING -> TRANSIENT (Duration: 30.00s)
11:40:08 - INFO - -------------------- 1% | 0.0s<0.1s | 3148.0 it/s
11:40:08 - INFO - ####---------------- 20% | 0.0s<0.1s | 3845.2 it/s
11:40:08 - INFO - ########------------ 40% | 0.0s<0.0s | 3794.2 it/s
11:40:08 - INFO - ############-------- 60% | 0.1s<0.0s | 3829.9 it/s
11:40:08 - INFO - ################---- 80% | 0.1s<0.0s | 3793.3 it/s
11:40:08 - INFO - #################### 100% | 0.1s<--:-- | 3790.1 it/s
11:40:08 - INFO - FINISHED -> TRANSIENT (total steps: 415, successful: 351, runtime: 114.38 ms)
[5]:
{'total_steps': 415, 'successful_steps': 351, 'runtime_ms': 114.37770500015176}
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)
11:40:08 - INFO - LOGGING (log: True)
11:40:08 - INFO - BLOCKS (total: 2, dynamic: 1, static: 1, eventful: 0)
11:40:08 - INFO - GRAPH (nodes: 2, edges: 2, alg. depth: 1, loop depth: 0, runtime: 0.084ms)
11:40:08 - INFO - STARTING -> TRANSIENT (Duration: 1500.00s)
11:40:08 - INFO - FINISHED -> TRANSIENT (total steps: 0, successful: 0, runtime: 2.80 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/stable/lib/python3.13/site-packages/pathsim/simulation.py:1744, in Simulation.run(self, duration, reset, adaptive)
1742 #enter tracker context and consume the run loop
1743 with tracker:
-> 1744 for _ in self._run_loop(duration, reset, adaptive, tracker=tracker):
1745 pass
1747 return tracker.stats
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/stable/lib/python3.13/site-packages/pathsim/simulation.py:1661, in Simulation._run_loop(self, duration, reset, adaptive, tracker)
1657 while self.time < end_time and self._active:
1658
1659 #advance the simulation by one (effective) timestep '_dt'
1660 try:
-> 1661 success, error_norm, scale, *_ = self.timestep(
1662 dt=_dt,
1663 adaptive=_adaptive
1664 )
1665 except StopSimulation as e:
1666 self.logger.info(f"STOP (StopSimulation: '{e}')")
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/stable/lib/python3.13/site-packages/pathsim/simulation.py:1476, in Simulation.timestep(self, dt, adaptive)
1472 for time_stage in self.engine.stages(self.time, dt):
1474 if is_implicit:
1475 #implicit: solve update equation (contains _update internally)
-> 1476 success, evals, solver_its = self._solve(time_stage, dt)
1477 total_evals += evals
1478 total_solver_its += solver_its
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/stable/lib/python3.13/site-packages/pathsim/simulation.py:1120, in Simulation._solve(self, t, dt)
1118 if not block:
1119 continue
-> 1120 self._solve_tracker.record(block, block.solve(t, dt))
1122 #check for convergence
1123 if self._solve_tracker.converged(self.tolerance_fpi):
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/stable/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/stable/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/stable/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/stable/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/stable/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/stable/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