Turning the knobs¶
A parameter study in a Modelica tool is usually a dialog box. Here the flat model is an
ordinary frozen dataclass, so a parameter study is a for loop over copies of it — and a
slider is six lines of ipywidgets.
This notebook is meant to be run. The figures below are drawn by matplotlib when the site
is built, so the page shows them; the sliders need a live kernel, so here they are a snapshot
of their last state. The file itself is attached to this page —
exploring.ipynb — but it reads its models from examples/, so it wants a
checkout rather than a bare download:
git clone https://gitlab.com/jorgeecardona/pymodelica && cd pymodelica
pip install "modelica[sim]" matplotlib ipywidgets
jupyter lab docs/notebooks/exploring.ipynb
from pathlib import Path
# The models come from `examples/` in the repository rather than from strings in here: a copy
# pasted into a notebook is a copy that goes stale, and the reason this page is executed when the
# site is built is that it then cannot show output from a model the project does not have.
root = next(
(p for p in [Path.cwd(), *Path.cwd().parents] if (p / "examples" / "models").is_dir()), None
)
if root is None:
message = (
"This notebook reads its models from `examples/` in the pymodelica repository. "
"Run it from a checkout:\n"
" git clone https://gitlab.com/jorgeecardona/pymodelica\n"
" cd pymodelica && jupyter lab docs/notebooks/"
)
raise RuntimeError(message)
One style block, used by every figure in this notebook. The colours are a validated categorical palette; the surface stays light in both site themes, because a figure is a picture and a picture does not get to change colour halfway down the page.
import matplotlib as mpl
import matplotlib.pyplot as plt
SURFACE, INK, MUTED = "#fcfcfb", "#0b0b0b", "#52514e"
BLUE, ORANGE = "#2a78d6", "#eb6834"
RAMP = ["#86b6ef", "#5598e7", "#2a78d6", "#1c5cab", "#104281"] # light -> dark, one hue
mpl.rcParams.update(
{
"figure.figsize": (8.4, 3.4),
"figure.dpi": 110,
"figure.facecolor": SURFACE,
"axes.facecolor": SURFACE,
"savefig.facecolor": SURFACE,
"axes.edgecolor": "#d4d3cd",
"axes.labelcolor": MUTED,
"axes.spines.top": False,
"axes.spines.right": False,
"axes.grid": True,
"axes.titlelocation": "left",
"axes.titlepad": 10,
"axes.titlesize": 11,
"grid.color": "#e8e7e2",
"grid.linewidth": 0.8,
"lines.linewidth": 2.0,
"legend.frameon": False,
"text.color": INK,
"xtick.color": MUTED,
"ytick.color": MUTED,
"font.size": 10,
}
)
%config InlineBackend.figure_formats = ["svg"]
The model, and what it has to turn¶
MassSpringDamper is the mechanical twin of the RLC circuit — a mass on a spring and a
damper, released from one metre. Flattening it resolves the connect equations and
qualifies every name, so the parameters come out as spring.c, damper.d, mass.m.
from modelica.lang import parse, unparse_expr
from modelica.sim import causalize, flatten, simulate
sources = [
parse((root / "examples" / "library" / "Translational.mo").read_text()),
parse((root / "examples" / "models" / "MassSpringDamper.mo").read_text()),
]
flat = flatten(sources, "MassSpringDamper")
for p in flat.parameters:
print(f"{p.name:16} = {unparse_expr(p.binding):<5} {p.unit:6} {p.comment}")
wall.s0 = 0 m Position it is fixed at spring.c = 4 N/m Stiffness spring.s_rel0 = 0 m Unstretched length damper.d = 0.4 N.s/m Damping constant mass.m = 1 kg Mass
Rewriting a binding¶
A parameter's value is the binding field of its FlatVariable, and a FlatVariable is
a frozen dataclass — so dataclasses.replace is the whole API. No re-parsing, no string
templating of .mo text, and the value that lands in the generated code is the one you
passed.
Causalization runs again on the copy. It has to: parameter values are folded into the emitted arithmetic rather than read from an array at run time, which is what makes the generated right-hand side straight-line code.
The cell below ends on a Result, and a result draws itself — so the chart that appears is
the library's own, not matplotlib's. There is more about it at the end of this notebook.
from dataclasses import replace
from modelica.ir import nodes as ir
from modelica.ir.flat import FlatModel
def tuned(model: FlatModel, values: dict[str, float]) -> FlatModel:
"""A copy of `model` with those parameters bound to those numbers."""
unknown = set(values) - {v.name for v in model.variables}
if unknown:
raise KeyError(f"{model.name} has no {sorted(unknown)}")
return replace(
model,
variables=tuple(
replace(v, binding=ir.RealLiteral(float(values[v.name]))) if v.name in values else v
for v in model.variables
),
)
def run(values: dict[str, float], stop: float = 12.0):
"""Tune, causalize, integrate — the three lines this notebook keeps repeating."""
return simulate(causalize(tuned(flat, values)), stop=stop, points=600)
run({"damper.d": 1.0})
A sweep is a loop¶
With mass.m = 1 and spring.c = 4, critical damping sits at
$d_c = 2\sqrt{cm} = 4$. Below it the mass overshoots and rings; above it, it crawls back
without crossing zero. One line per damping value, dark for more damping:
DAMPING = [0.4, 1.0, 2.0, 4.0, 8.0]
fig, ax = plt.subplots()
for d, colour in zip(DAMPING, RAMP, strict=True):
result = run({"damper.d": d})
ax.plot(result.time, result["mass.s"], color=colour, label=f"d = {d}")
ax.axhline(0, color="#d4d3cd", linewidth=1, zorder=0)
ax.set_title("Displacement of the mass, damping 0.4 → 8 N·s/m (critical damping is 4)")
ax.set_xlabel("time [s]")
ax.set_ylabel("mass.s [m]")
ax.legend(loc="center left", bbox_to_anchor=(1.02, 0.5), labelcolor=MUTED)
fig.subplots_adjust(right=0.84)
plt.show()
The d = 4 curve is the fastest one that never goes below zero — which is the definition
of critical damping, recovered from the simulator rather than assumed by it.
Sweeping finer turns the same picture into two numbers per run: how far the mass overshoots, and how long it takes to stay inside ±2 cm.
import numpy as np
BAND = 0.02
sweep = np.linspace(0.2, 10.0, 30)
overshoot, settling = [], []
for d in sweep:
result = run({"damper.d": float(d)})
s = result["mass.s"]
overshoot.append(max(0.0, -s.min()))
# The last sample outside the band, then interpolated to the crossing itself -- reading
# the sample time instead would quantize this curve to the output grid.
outside = np.flatnonzero(np.abs(s) > BAND)
last = int(outside[-1]) if len(outside) else 0
here, then = abs(float(s[last])), abs(float(s[min(last + 1, len(s) - 1)]))
span = float(result.time[min(last + 1, len(s) - 1)] - result.time[last])
settling.append(
float(result.time[last]) + (span * (here - BAND) / (here - then) if here > then else 0.0)
)
fig, (left, right) = plt.subplots(1, 2, figsize=(8.4, 3.2))
left.plot(sweep, np.array(overshoot) * 100, color=BLUE)
left.set_title("Overshoot below rest")
left.set_ylabel("overshoot [cm]")
right.plot(sweep, settling, color=BLUE)
right.set_title("Time to stay inside ±2 cm")
right.set_ylabel("settling time [s]")
for ax in (left, right):
ax.axvline(4.0, color=ORANGE, linewidth=1.5, linestyle=(0, (4, 3)))
ax.annotate(
"critical, d = 4",
xy=(4.0, ax.get_ylim()[1]),
xytext=(4.5, ax.get_ylim()[1] * 0.92),
color=ORANGE,
fontsize=9,
)
ax.set_xlabel("damping d [N·s/m]")
plt.tight_layout()
plt.show()
Overshoot reaches zero at $d = 4$ and stays there. Settling time has its minimum just short of it — slightly underdamped settles fastest, which is why controllers are tuned there and not at the textbook value.
The steps in the settling curve are not solver noise. Settling time is the last moment the
trajectory is outside the band, so it jumps whenever the ringing dies down enough that an
earlier lobe becomes the last one to leave — a discrete change in which crossing is being
measured, in a quantity that is otherwise continuous in d.
The same thing, with sliders¶
run takes about eighty milliseconds for this model, so a slider can call it directly.
continuous_update=False means the resimulation happens when the handle is dropped
rather than on every pixel of the drag.
def draw(m: float, c: float, d: float) -> None:
"""Displacement and phase portrait for one (m, c, d), plus the damping ratio."""
result = run({"mass.m": m, "spring.c": c, "damper.d": d})
zeta = d / (2.0 * np.sqrt(c * m))
regime = "underdamped" if zeta < 1 else ("critically damped" if zeta == 1 else "overdamped")
fig, (left, right) = plt.subplots(1, 2, figsize=(8.4, 3.2))
left.plot(result.time, result["mass.s"], color=BLUE)
left.axhline(0, color="#d4d3cd", linewidth=1, zorder=0)
left.set_title(f"ζ = {zeta:.2f} — {regime}")
left.set_xlabel("time [s]")
left.set_ylabel("mass.s [m]")
left.set_ylim(-1.1, 1.1)
right.plot(result["mass.s"], result["mass.v"], color=BLUE)
right.plot([1.0], [0.0], marker="o", markersize=8, color=ORANGE)
right.annotate("release", xy=(1.0, 0.0), xytext=(0.55, 0.45), color=ORANGE, fontsize=9)
right.set_title("Phase portrait")
right.set_xlabel("mass.s [m]")
right.set_ylabel("mass.v [m/s]")
right.set_xlim(-1.2, 1.2)
right.set_ylim(-2.2, 2.2)
plt.tight_layout()
plt.show()
draw(m=1.0, c=4.0, d=0.4)
import ipywidgets as widgets
_ = widgets.interact(
draw,
m=widgets.FloatSlider(
value=1.0, min=0.2, max=4.0, step=0.1, description="mass.m", continuous_update=False
),
c=widgets.FloatSlider(
value=4.0, min=0.5, max=20.0, step=0.5, description="spring.c", continuous_update=False
),
d=widgets.FloatSlider(
value=0.4, min=0.0, max=10.0, step=0.2, description="damper.d", continuous_update=False
),
)
The spiral tightens as damper.d grows and opens into a single curve into the origin once
ζ passes 1. Set damper.d to exactly 0 and the phase portrait closes into an ellipse: no
dissipation, so the orbit repeats forever — and the ellipse not drifting outward over
twelve seconds is a statement about the integrator, not about the model.
Every variable, without matplotlib¶
Result.plot() is the library's own chart: SVG drawn by a few hundred lines of vanilla
JavaScript inside a sandboxed <iframe>, with no plotting dependency and no CDN behind
it. It stays interactive on this page — drag to zoom, click the legend, read values at the
cursor — because it does not need a kernel, only a browser.
Every variable of the flat model is in the legend, aliases included: spring.f,
damper.flange_b.s, the lot. The moving ones are ticked on to begin with.
run({"damper.d": 0.4}).plot()
Next: the shape of the equations — the same models as matrices, and what matching and BLT sorting do to them.