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