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_5314/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)
10:52:26 - INFO - LOGGING (log: True)
10:52:26 - INFO - BLOCKS (total: 2, dynamic: 1, static: 1, eventful: 0)
10:52:26 - INFO - GRAPH (nodes: 2, edges: 2, alg. depth: 1, loop depth: 0, runtime: 0.161ms)
10:52:26 - INFO - STARTING -> TRANSIENT (Duration: 30.00s)
10:52:26 - INFO - -------------------- 1% | 0.0s<0.2s | 1784.6 it/s
10:52:26 - INFO - ####---------------- 20% | 0.0s<0.1s | 2124.8 it/s
10:52:27 - INFO - ########------------ 40% | 0.1s<0.1s | 2169.4 it/s
10:52:27 - INFO - ############-------- 60% | 0.1s<0.1s | 2171.4 it/s
10:52:27 - INFO - ################---- 80% | 0.2s<0.0s | 2172.6 it/s
10:52:27 - INFO - #################### 100% | 0.2s<--:-- | 2087.8 it/s
10:52:27 - INFO - FINISHED -> TRANSIENT (total steps: 415, successful: 351, runtime: 207.78 ms)
[5]:
{'total_steps': 415, 'successful_steps': 351, 'runtime_ms': 207.78222900025867}
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)
10:52:27 - INFO - LOGGING (log: True)
10:52:27 - INFO - BLOCKS (total: 2, dynamic: 1, static: 1, eventful: 0)
10:52:27 - INFO - GRAPH (nodes: 2, edges: 2, alg. depth: 1, loop depth: 0, runtime: 0.130ms)
10:52:27 - INFO - STARTING -> TRANSIENT (Duration: 1500.00s)
10:52:27 - INFO - FINISHED -> TRANSIENT (total steps: 0, successful: 0, runtime: 4.89 ms)
---------------------------------------------------------------------------
ValueError 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.23.0/lib/python3.13/site-packages/pathsim/simulation.py:2032, in Simulation.run(self, duration, reset, adaptive)
2030 #enter tracker context and consume the run loop
2031 with tracker:
-> 2032 for _ in self._run_loop(duration, reset, adaptive, tracker=tracker):
2033 pass
2035 return tracker.stats
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.23.0/lib/python3.13/site-packages/pathsim/simulation.py:1949, in Simulation._run_loop(self, duration, reset, adaptive, tracker)
1945 while self.time < end_time and self._active:
1946
1947 #advance the simulation by one (effective) timestep '_dt'
1948 try:
-> 1949 success, error_norm, scale, *_ = self.timestep(
1950 dt=_dt,
1951 adaptive=_adaptive
1952 )
1953 except StopSimulation as e:
1954 self.logger.info(f"STOP (StopSimulation: '{e}')")
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.23.0/lib/python3.13/site-packages/pathsim/simulation.py:1759, in Simulation.timestep(self, dt, adaptive)
1755 for time_stage in self.engine.stages(self.time, dt):
1757 if is_implicit:
1758 #implicit: solve update equation (contains _update internally)
-> 1759 success, evals, solver_its = self._solve(time_stage, dt)
1760 total_evals += evals
1761 total_solver_its += solver_its
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.23.0/lib/python3.13/site-packages/pathsim/simulation.py:1124, in Simulation._solve(self, t, dt)
1122 if not block:
1123 continue
-> 1124 self._solve_tracker.record(block, block.solve(t, dt))
1126 #check for convergence
1127 if self._solve_tracker.converged(self.tolerance_fpi):
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.23.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.23.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.23.0/lib/python3.13/site-packages/pathsim/optim/anderson.py:394, in NewtonAnderson.step(self, x, g, jac)
391 return super().step(x, g)
392 else:
393 #newton step with residual
--> 394 _x, res_norm = self._newton(x, g, jac)
396 #anderson step with no residual
397 y, _ = super().step(_x, g)
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.23.0/lib/python3.13/site-packages/pathsim/optim/anderson.py:361, in NewtonAnderson._newton(self, x, g, jac)
358 _lu = lu_factor(_A)
359 self._A, self._lu = _A.copy(), _lu
--> 361 return _x - lu_solve(_lu, _res), np.linalg.norm(_res)
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.23.0/lib/python3.13/site-packages/scipy/linalg/_decomp_lu.py:188, in lu_solve(lu_and_piv, b, trans, overwrite_b, check_finite)
135 """Solve an equation system, a x = b, given the LU factorization of a
136
137 The documentation is written assuming array arguments are of specified
(...) 185
186 """
187 (lu, piv) = lu_and_piv
--> 188 return _lu_solve(lu, piv, b, trans=trans, overwrite_b=overwrite_b,
189 check_finite=check_finite)
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.23.0/lib/python3.13/site-packages/scipy/_lib/_util.py:1181, in _apply_over_batch.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
1179 # Early exit if call is not batched
1180 if not any(batch_shapes):
-> 1181 return f(*arrays, *other_args, **kwargs)
1183 # Determine broadcasted batch shape
1184 batch_shape = np.broadcast_shapes(*batch_shapes) # Gives OK error message
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.23.0/lib/python3.13/site-packages/scipy/linalg/_decomp_lu.py:195, in _lu_solve(lu, piv, b, trans, overwrite_b, check_finite)
192 @_apply_over_batch(('lu', 2), ('piv', 1), ('b', '1|2'))
193 def _lu_solve(lu, piv, b, trans, overwrite_b, check_finite):
194 if check_finite:
--> 195 b1 = asarray_chkfinite(b)
196 else:
197 b1 = asarray(b)
File ~/checkouts/readthedocs.org/user_builds/pathsim/envs/v0.23.0/lib/python3.13/site-packages/numpy/lib/function_base.py:630, in asarray_chkfinite(a, dtype, order)
628 a = asarray(a, dtype=dtype, order=order)
629 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
--> 630 raise ValueError(
631 "array must not contain infs or NaNs")
632 return a
ValueError: array must not contain infs or NaNs
[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