{
"cells": [
{
"cell_type": "markdown",
"id": "73def54e",
"metadata": {},
"source": [
"# Roasting a batch\n",
"\n",
"Every other notebook here simulates something the project made up. This one simulates a\n",
"machine that exists: a **Loring S15**, fifteen kilograms of green coffee at a time, which\n",
"roasted five batches on the evening of 30 July 2026 while a bridge logged it over Modbus\n",
"every one and a half seconds.\n",
"\n",
"So there is a right answer, and it was not computed — it was measured. The figures below\n",
"put the model's curve on top of the roaster's own bean probe and leave the gap visible.\n",
"\n",
"One of the four states is the thermocouple itself, and that matters more than it sounds. A\n",
"bean probe is a lump of steel with its own heat capacity, so for the first minute of a\n",
"roast it is not measuring the coffee so much as catching up with it — which is the whole\n",
"explanation of the *turning point*, the dive every roaster watches for in the first\n",
"seconds.\n",
"\n",
"This notebook is meant to be **run**. The figures are drawn by matplotlib when the site is\n",
"built; the sliders at the end need a live kernel, so on this page they are a snapshot of\n",
"their last state. The file is attached — [roasting.ipynb](roasting.ipynb) — but it reads\n",
"the model from `examples/`, so it wants a 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/roasting.ipynb\n",
"```\n",
"\n",
"The measurement it checks against is carried inline, a little further down. The roast logs\n",
"live in a different repository and this page does not want to depend on one."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "220ed581",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\"Seconds after the charge beyond which the probe's arrival transient is over.\\n\\nBefore it the thermocouple is still catching up with the bean mass that buried it, and on\\na batch charged straight off a drop its reading at the mark is contaminated. The model\\ncovers that stretch -- it is what `T_probe` is for -- but it is where the disagreement\\nlives, so the figures report either side of it separately.\\n\""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from pathlib import Path\n",
"\n",
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"\n",
"# The model comes from `examples/` in the repository rather than from a string 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 model 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)\n",
"\n",
"SURFACE, INK, MUTED, LINE = \"#fcfcfb\", \"#0b0b0b\", \"#52514e\", \"#d4d3cd\"\n",
"BLUE, ORANGE, AQUA, YELLOW = \"#2a78d6\", \"#eb6834\", \"#1baf7a\", \"#eda100\"\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\": LINE,\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\"]\n",
"\n",
"\n",
"def minutes(ax):\n",
" \"\"\"Roasts are read in minutes and seconds, not in seconds.\"\"\"\n",
" ax.set_xlabel(\"time from charge (min)\")\n",
" ax.xaxis.set_major_formatter(lambda s, _: f\"{int(s // 60)}:{int(s % 60):02d}\")\n",
" ax.xaxis.set_major_locator(mpl.ticker.MultipleLocator(120))\n",
"\n",
"\n",
"def label(ax, x, y, text, colour, dy=0):\n",
" \"\"\"Direct labels, one per series.\n",
"\n",
" Four series here sit under 3:1 against the surface, so every one of them is named on the\n",
" plot rather than left to colour alone.\n",
" \"\"\"\n",
" ax.annotate(\n",
" text,\n",
" xy=(x, y),\n",
" xytext=(6, dy),\n",
" textcoords=\"offset points\",\n",
" color=colour,\n",
" fontsize=9,\n",
" va=\"center\",\n",
" )\n",
"\n",
"\n",
"SETTLED = 120.0\n",
"\"\"\"Seconds after the charge beyond which the probe's arrival transient is over.\n",
"\n",
"Before it the thermocouple is still catching up with the bean mass that buried it, and on\n",
"a batch charged straight off a drop its reading at the mark is contaminated. The model\n",
"covers that stretch -- it is what `T_probe` is for -- but it is where the disagreement\n",
"lives, so the figures report either side of it separately.\n",
"\"\"\""
]
},
{
"cell_type": "markdown",
"id": "bb07f45b",
"metadata": {},
"source": [
"## The machine\n",
"\n",
"Four states, and each is a different kind of thing:\n",
"\n",
"| state | what it is |\n",
"|---|---|\n",
"| `T_mach` | everything hot that is not coffee — steel, lining, recirculated air |\n",
"| `T_bean` | the coffee |\n",
"| `T_probe` | the thermocouple *in* the coffee, which is not the same as the coffee |\n",
"| `m_w` | the water still in the beans |\n",
"\n",
"The burner pushes heat into the machine, the machine leaks some to the room and hands the\n",
"rest to the coffee, the coffee boils water off once it is past 100 °C, and the probe trails\n",
"along behind.\n",
"\n",
"`moisture` is a **parameter**, not a fitted number, and it does real work: the charge's heat\n",
"capacity is `m_dry*cp_dry + m_w*cp_water`, so it falls as the batch dries. `m_dry` and the\n",
"starting water are computed *from* it — a parameter bound to an expression over other\n",
"parameters, resolved to a number before any integration happens.\n",
"\n",
"The nine-level chain on `burner` is not invented. It is this batch's own recorded gas\n",
"trace, reduced to the steps worth carrying."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "750226b5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"model RoastingMachine \"A Loring S15 coffee roaster over one batch -- gas on a clock, water on a threshold\"\n",
" parameter Real T_amb(unit = \"K\") = 293.15 \"The roastery, 20 degC\";\n",
" parameter Real m_green(unit = \"kg\") = 15 \"The charge. Rated, not weighed\";\n",
" parameter Real moisture = 0.115 \"Water in the green coffee, by mass\";\n",
" parameter Real m_dry(unit = \"kg\") = m_green * (1 - moisture) \"Dry matter\";\n",
" parameter Real cp_dry(unit = \"J/(kg.K)\") = 1800 \"Dry coffee\";\n",
" parameter Real cp_water(unit = \"J/(kg.K)\") = 4186 \"Water\";\n",
" parameter Real L(unit = \"J/kg\") = 2260000.0 \"Latent heat of vaporisation\";\n",
" parameter Real k_dry(unit = \"1/(K.s)\") = 4.67e-6 \"Drying rate above boiling\";\n",
" parameter Real C_mach(unit = \"J/K\") = 530000.0 \"Machine hot mass: steel, lining, loop air\";\n",
" parameter Real UA(unit = \"W/K\") = 261 \"Air-to-bean coupling over the whole bed\";\n",
" parameter Real G_loss(unit = \"W/K\") = 259 \"Loop to ambient: the bleed and the shell\";\n",
" parameter Real Q_max(unit = \"W\") = 139000.0 \"Burner at 100 percent\";\n",
" parameter Real tau_probe(unit = \"s\") = 17.8 \"Bean probe time constant\";\n",
" Real T_mach(unit = \"K\", start = 428.23, fixed = true) \"Machine temperature\";\n",
" Real T_bean(unit = \"K\", start = 293.15, fixed = true) \"Bean temperature -- the coffee itself\";\n",
" Real T_probe(unit = \"K\", start = 489.98, fixed = true) \"Bean probe temperature\";\n",
" Real m_w(unit = \"kg\", start = m_green * moisture, fixed = true) \"Water still in the beans\";\n",
" Real C_bean(unit = \"J/K\") \"Heat capacity of the charge -- falls as it dries\";\n",
" Real evap(unit = \"kg/s\") \"Evaporation rate\";\n",
" Real burner \"Gas modulation, per cent -- the machine's own signal\";\n",
" Real Q_burner(unit = \"W\") \"Burner heat into the loop\";\n",
" Real bean(unit = \"degC\") \"What the coffee is at, in degC\";\n",
" Real probe(unit = \"degC\") \"What the roaster reads, in degC\";\n",
" Real lag(unit = \"K\") \"How far the coffee runs ahead of its own probe\";\n",
" Real loss \"Weight lost as steam, per cent of the charge\";\n",
"equation\n",
" burner = if time < 66.7 then 20 elseif time < 93.0 then 27 elseif time < 526.9 then 85 elseif time < 618.0 then 79 elseif time < 688.5 then 70 elseif time < 714.6 then 61 elseif time < 754.4 then 50 elseif time < 780.8 then 37 else 20;\n",
" Q_burner = burner / 100 * Q_max;\n",
" C_mach * der(T_mach) = Q_burner - UA * (T_mach - T_bean) - G_loss * (T_mach - T_amb);\n",
" evap = if T_bean > 373.15 then k_dry * m_w * (T_bean - 373.15) else 0;\n",
" C_bean = m_dry * cp_dry + m_w * cp_water;\n",
" C_bean * der(T_bean) = UA * (T_mach - T_bean) - L * evap;\n",
" der(m_w) = -evap;\n",
" tau_probe * der(T_probe) = T_bean - T_probe;\n",
" bean = T_bean - 273.15;\n",
" probe = T_probe - 273.15;\n",
" lag = T_bean - T_probe;\n",
" loss = 100 * (m_green * moisture - m_w) / m_green;\n",
"end RoastingMachine;\n",
"\n"
]
}
],
"source": [
"print((root / \"examples\" / \"models\" / \"RoastingMachine.mo\").read_text())"
]
},
{
"cell_type": "markdown",
"id": "8cde9f71",
"metadata": {},
"source": [
"## One batch\n",
"\n",
"`parse` → `flatten` → `causalize` → `simulate`, the same four calls as every other model in\n",
"the corpus. What comes back is fourteen and a half minutes of roast."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "6143846a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"states ['T_bean', 'T_mach', 'T_probe', 'm_w']\n",
"events 9\n",
"drop the probe reads 213.5 degC; the coffee is at 213.7 degC\n",
"water 1.324 kg left, 2.67 % of the charge driven off as steam\n"
]
}
],
"source": [
"from modelica.lang import parse\n",
"from modelica.sim import causalize, flatten, simulate\n",
"\n",
"flat = flatten(\n",
" [parse((root / \"examples\" / \"models\" / \"RoastingMachine.mo\").read_text())], \"RoastingMachine\"\n",
")\n",
"system = causalize(flat)\n",
"\n",
"STOP = 863.6\n",
"result = simulate(system, stop=STOP, tolerance=1e-10, points=900)\n",
"\n",
"end = result.at(STOP)\n",
"print(f\"states {sorted(system.states)}\")\n",
"print(f\"events {len(result.events)}\")\n",
"print(f\"drop the probe reads {end['probe']:.1f} degC; the coffee is at {end['bean']:.1f} degC\")\n",
"print(f\"water {end['m_w']:.3f} kg left, {end['loss']:.2f} % of the charge driven off as steam\")"
]
},
{
"cell_type": "markdown",
"id": "cf9b13e0",
"metadata": {},
"source": [
"## Two kinds of event in one model\n",
"\n",
"Nine events, and they are not the same species.\n",
"\n",
"**Eight are on the clock.** A gas step fires because it is twenty past ten — the condition\n",
"is on `time` itself, so every one of those instants is known before the integration starts.\n",
"\n",
"**One is on the state.** `evap` is zero until the beans pass 100 °C, and *when* that happens\n",
"depends on the trajectory. Nobody knows it in advance; the solver has to hunt for it, the\n",
"way `DrainingTank` hunts for the moment a tank runs dry. It lands at 111.4 s.\n",
"\n",
"Having both kinds in one model is why this example is in the corpus. The burner gets its own\n",
"panel rather than a second y-axis — per cent and degrees are different measures and do not\n",
"belong on one scale."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "2444cd2d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" t = 66.7 s ( 1.11 min) time < 66.7 changed\n",
" t = 93.0 s ( 1.55 min) time < 93.0 changed\n",
" t = 111.4 s ( 1.86 min) T_bean > 373.15 changed\n",
" t = 526.9 s ( 8.78 min) time < 526.9 changed\n",
" t = 618.0 s (10.30 min) time < 618.0 changed\n",
" t = 688.5 s (11.47 min) time < 688.5 changed\n",
" t = 714.6 s (11.91 min) time < 714.6 changed\n",
" t = 754.4 s (12.57 min) time < 754.4 changed\n",
" t = 780.8 s (13.01 min) time < 780.8 changed\n"
]
}
],
"source": [
"for event in result.events:\n",
" print(f\" t = {event.time:6.1f} s ({event.time / 60:5.2f} min) {event.cause}\")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "9a4b7fc3",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"t = result.time\n",
"bean = result.values[result.names.index(\"bean\")]\n",
"probe = result.values[result.names.index(\"probe\")]\n",
"mach = result.values[result.names.index(\"T_mach\")] - 273.15\n",
"gas = result.values[result.names.index(\"burner\")]\n",
"water = result.values[result.names.index(\"m_w\")]\n",
"\n",
"fig, (top, bottom) = plt.subplots(\n",
" 2, 1, figsize=(8.4, 5.0), sharex=True, height_ratios=[2.4, 1], gridspec_kw={\"hspace\": 0.15}\n",
")\n",
"\n",
"top.plot(t, mach, color=AQUA)\n",
"top.plot(t, bean, color=BLUE)\n",
"top.plot(t, probe, color=ORANGE, linewidth=1.5, linestyle=(0, (5, 2)))\n",
"label(top, t[-1], mach[-1], \"machine\", AQUA, dy=9)\n",
"label(top, t[-1], bean[-1], \"coffee\", BLUE, dy=-4)\n",
"label(top, t[-1], probe[-1], \"probe\", ORANGE, dy=-16)\n",
"top.set_title(\"Machine, coffee, and the probe that reports on the coffee\")\n",
"top.set_ylabel(\"degC\")\n",
"top.set_xlim(0, STOP * 1.09)\n",
"\n",
"bottom.plot(t, gas, color=YELLOW, drawstyle=\"steps-post\")\n",
"label(bottom, t[-1], gas[-1], \"burner\", YELLOW)\n",
"bottom.set_ylabel(\"%\")\n",
"bottom.set_ylim(0, 100)\n",
"minutes(bottom)\n",
"\n",
"boil = [e.time for e in result.events if \"T_bean\" in e.cause] # the one state event\n",
"for event in result.events:\n",
" on_state = event.time in boil\n",
" for ax in (top, bottom):\n",
" ax.axvline(\n",
" event.time,\n",
" color=BLUE if on_state else LINE,\n",
" linewidth=1.2 if on_state else 0.8,\n",
" linestyle=\":\" if on_state else \"-\",\n",
" zorder=0,\n",
" )\n",
"top.annotate(\n",
" \"beans pass the drying onset:\\nthe water starts leaving\",\n",
" xy=(boil[0], 40.0),\n",
" xytext=(215, 38),\n",
" fontsize=8.5,\n",
" color=MUTED,\n",
" arrowprops={\"arrowstyle\": \"-\", \"color\": MUTED, \"linewidth\": 0.8},\n",
")\n",
"\n",
"fig.legend(\n",
" handles=[\n",
" mpl.lines.Line2D([], [], color=AQUA, label=\"machine\"),\n",
" mpl.lines.Line2D([], [], color=BLUE, label=\"coffee\"),\n",
" mpl.lines.Line2D([], [], color=ORANGE, linestyle=(0, (5, 2)), label=\"bean probe\"),\n",
" mpl.lines.Line2D([], [], color=YELLOW, label=\"burner\"),\n",
" mpl.lines.Line2D([], [], color=LINE, linewidth=0.8, label=\"time event\"),\n",
" mpl.lines.Line2D([], [], color=BLUE, linewidth=1.2, linestyle=\":\", label=\"state event\"),\n",
" ],\n",
" loc=\"lower center\",\n",
" ncol=6,\n",
" bbox_to_anchor=(0.5, -0.05),\n",
" fontsize=9,\n",
")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "62acd55a",
"metadata": {},
"source": [
"## What the roaster actually read\n",
"\n",
"Here is the measurement. Five batches, downsampled to fifteen-second spacing and carried\n",
"inline so this page owes nothing to another repository: for each one, the gas trace the\n",
"machine followed and what its bean probe reported."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "bade1bbc",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"5 batches, 290 points\n",
" batch 1: 14.39 min, 9 gas steps, turn 46.8 s, drop 216.06 degC\n",
" batch 2: 13.90 min, 28 gas steps, turn 43.8 s, drop 216.94 degC\n",
" batch 3: 14.27 min, 24 gas steps, turn 44.3 s, drop 215.72 degC\n",
" batch 4: 14.18 min, 21 gas steps, turn 42.3 s, drop 215.89 degC\n",
" batch 5: 13.98 min, 10 gas steps, turn 41.0 s, drop 215.11 degC\n"
]
}
],
"source": [
"BATCHES = [\n",
" # batch 1: 14.39 min, 9 gas steps\n",
" {\n",
" \"mach0\": 237.38,\n",
" \"stop\": 863.63,\n",
" \"turn\": 46.75,\n",
" \"gas\": [\n",
" (0.0, 20.0),\n",
" (66.7, 27.01),\n",
" (93.0, 84.97),\n",
" (526.9, 79.29),\n",
" (618.0, 70.0),\n",
" (688.5, 61.46),\n",
" (714.6, 50.0),\n",
" (754.4, 37.42),\n",
" (780.8, 20.0),\n",
" ],\n",
" \"bt\": [\n",
" (0.0, 216.83),\n",
" (15.0, 102.84),\n",
" (30.0, 70.45),\n",
" (45.0, 64.49),\n",
" (60.0, 67.25),\n",
" (75.0, 72.73),\n",
" (90.0, 79.94),\n",
" (105.0, 87.57),\n",
" (120.0, 94.87),\n",
" (135.0, 101.41),\n",
" (150.0, 107.61),\n",
" (165.0, 113.39),\n",
" (180.0, 118.96),\n",
" (195.0, 123.84),\n",
" (210.0, 127.71),\n",
" (225.0, 131.7),\n",
" (240.0, 135.39),\n",
" (255.0, 138.73),\n",
" (270.0, 141.76),\n",
" (285.0, 144.94),\n",
" (300.0, 147.65),\n",
" (315.0, 150.33),\n",
" (330.0, 153.13),\n",
" (345.0, 155.73),\n",
" (360.0, 158.27),\n",
" (375.0, 161.0),\n",
" (390.0, 163.16),\n",
" (405.0, 165.6),\n",
" (420.0, 168.11),\n",
" (435.0, 170.18),\n",
" (450.0, 172.4),\n",
" (465.0, 174.37),\n",
" (480.0, 176.53),\n",
" (495.0, 178.58),\n",
" (510.0, 180.56),\n",
" (525.0, 182.76),\n",
" (540.0, 184.74),\n",
" (555.0, 186.74),\n",
" (570.0, 188.74),\n",
" (585.0, 190.57),\n",
" (600.0, 192.48),\n",
" (615.0, 194.79),\n",
" (630.0, 196.86),\n",
" (645.0, 198.89),\n",
" (660.0, 201.03),\n",
" (675.0, 203.12),\n",
" (690.0, 205.32),\n",
" (705.0, 207.26),\n",
" (720.0, 208.83),\n",
" (735.0, 210.03),\n",
" (750.0, 211.07),\n",
" (765.0, 211.94),\n",
" (780.0, 212.28),\n",
" (795.0, 212.68),\n",
" (810.0, 213.5),\n",
" (825.0, 214.15),\n",
" (840.0, 214.68),\n",
" (855.0, 215.56),\n",
" (863.6, 216.06),\n",
" ],\n",
" },\n",
" # batch 2: 13.90 min, 28 gas steps\n",
" {\n",
" \"mach0\": 234.67,\n",
" \"stop\": 833.7,\n",
" \"turn\": 43.83,\n",
" \"gas\": [\n",
" (0.0, 20.0),\n",
" (50.0, 74.5),\n",
" (76.2, 27.19),\n",
" (102.2, 50.99),\n",
" (128.1, 100.0),\n",
" (154.2, 78.47),\n",
" (190.9, 84.56),\n",
" (227.6, 90.58),\n",
" (256.8, 98.53),\n",
" (283.0, 92.51),\n",
" (308.8, 100.0),\n",
" (334.7, 93.85),\n",
" (360.8, 99.59),\n",
" (387.0, 93.21),\n",
" (413.1, 100.0),\n",
" (486.5, 94.91),\n",
" (512.8, 100.0),\n",
" (538.9, 91.92),\n",
" (565.2, 78.18),\n",
" (591.3, 64.21),\n",
" (623.3, 54.44),\n",
" (654.0, 49.06),\n",
" (680.3, 38.88),\n",
" (706.3, 67.48),\n",
" (732.5, 55.43),\n",
" (759.5, 50.11),\n",
" (785.2, 37.54),\n",
" (811.0, 20.0),\n",
" ],\n",
" \"bt\": [\n",
" (0.0, 192.17),\n",
" (15.0, 100.61),\n",
" (30.0, 68.4),\n",
" (45.0, 62.76),\n",
" (60.0, 66.1),\n",
" (75.0, 73.07),\n",
" (90.0, 80.12),\n",
" (105.0, 86.27),\n",
" (120.0, 92.48),\n",
" (135.0, 99.1),\n",
" (150.0, 105.64),\n",
" (165.0, 111.43),\n",
" (180.0, 116.5),\n",
" (195.0, 121.2),\n",
" (210.0, 125.92),\n",
" (225.0, 129.94),\n",
" (240.0, 133.96),\n",
" (255.0, 137.77),\n",
" (270.0, 141.22),\n",
" (285.0, 144.62),\n",
" (300.0, 147.44),\n",
" (315.0, 150.52),\n",
" (330.0, 153.85),\n",
" (345.0, 156.5),\n",
" (360.0, 159.14),\n",
" (375.0, 161.82),\n",
" (390.0, 164.49),\n",
" (405.0, 167.17),\n",
" (420.0, 169.58),\n",
" (435.0, 172.13),\n",
" (450.0, 174.57),\n",
" (465.0, 177.0),\n",
" (480.0, 179.52),\n",
" (495.0, 181.95),\n",
" (510.0, 184.26),\n",
" (525.0, 186.59),\n",
" (540.0, 189.03),\n",
" (555.0, 191.16),\n",
" (570.0, 193.3),\n",
" (585.0, 195.54),\n",
" (600.0, 197.32),\n",
" (615.0, 198.64),\n",
" (630.0, 200.58),\n",
" (645.0, 202.18),\n",
" (660.0, 203.49),\n",
" (675.0, 204.98),\n",
" (690.0, 206.33),\n",
" (705.0, 207.51),\n",
" (720.0, 209.04),\n",
" (735.0, 210.66),\n",
" (750.0, 211.87),\n",
" (765.0, 212.57),\n",
" (780.0, 213.68),\n",
" (795.0, 214.65),\n",
" (810.0, 215.5),\n",
" (825.0, 216.35),\n",
" (833.7, 216.94),\n",
" ],\n",
" },\n",
" # batch 3: 14.27 min, 24 gas steps\n",
" {\n",
" \"mach0\": 234.23,\n",
" \"stop\": 856.39,\n",
" \"turn\": 44.26,\n",
" \"gas\": [\n",
" (0.0, 20.0),\n",
" (50.3, 87.48),\n",
" (75.8, 20.0),\n",
" (101.3, 30.99),\n",
" (127.0, 100.0),\n",
" (152.6, 92.63),\n",
" (178.2, 79.76),\n",
" (205.3, 69.06),\n",
" (230.9, 100.0),\n",
" (273.0, 92.63),\n",
" (301.7, 99.23),\n",
" (429.5, 91.52),\n",
" (455.1, 100.0),\n",
" (504.6, 93.56),\n",
" (530.1, 99.29),\n",
" (555.8, 74.15),\n",
" (581.3, 79.76),\n",
" (606.9, 53.97),\n",
" (632.6, 44.26),\n",
" (658.3, 51.46),\n",
" (683.7, 43.97),\n",
" (709.4, 60.76),\n",
" (735.1, 34.97),\n",
" (786.8, 22.45),\n",
" ],\n",
" \"bt\": [\n",
" (0.0, 171.56),\n",
" (15.0, 90.21),\n",
" (30.0, 64.74),\n",
" (45.0, 61.55),\n",
" (60.0, 65.8),\n",
" (75.0, 73.38),\n",
" (90.0, 80.39),\n",
" (105.0, 86.41),\n",
" (120.0, 92.25),\n",
" (135.0, 98.54),\n",
" (150.0, 105.35),\n",
" (165.0, 111.05),\n",
" (180.0, 116.35),\n",
" (195.0, 121.25),\n",
" (210.0, 125.76),\n",
" (225.0, 129.51),\n",
" (240.0, 133.19),\n",
" (255.0, 137.1),\n",
" (270.0, 140.84),\n",
" (285.0, 144.41),\n",
" (300.0, 147.29),\n",
" (315.0, 150.17),\n",
" (330.0, 153.32),\n",
" (345.0, 156.22),\n",
" (360.0, 158.83),\n",
" (375.0, 161.53),\n",
" (390.0, 164.25),\n",
" (405.0, 166.9),\n",
" (420.0, 169.52),\n",
" (435.0, 172.23),\n",
" (450.0, 174.46),\n",
" (465.0, 176.73),\n",
" (480.0, 179.22),\n",
" (495.0, 181.66),\n",
" (510.0, 184.04),\n",
" (525.0, 186.28),\n",
" (540.0, 188.73),\n",
" (555.0, 191.02),\n",
" (570.0, 192.94),\n",
" (585.0, 195.19),\n",
" (600.0, 196.97),\n",
" (615.0, 198.83),\n",
" (630.0, 200.38),\n",
" (645.0, 201.69),\n",
" (660.0, 203.4),\n",
" (675.0, 204.85),\n",
" (690.0, 206.14),\n",
" (705.0, 207.46),\n",
" (720.0, 208.56),\n",
" (735.0, 209.44),\n",
" (750.0, 210.11),\n",
" (765.0, 210.82),\n",
" (780.0, 211.61),\n",
" (795.0, 212.41),\n",
" (810.0, 212.96),\n",
" (825.0, 213.83),\n",
" (840.0, 214.58),\n",
" (855.0, 215.47),\n",
" (856.4, 215.72),\n",
" ],\n",
" },\n",
" # batch 4: 14.18 min, 21 gas steps\n",
" {\n",
" \"mach0\": 233.63,\n",
" \"stop\": 851.0,\n",
" \"turn\": 42.26,\n",
" \"gas\": [\n",
" (0.0, 20.0),\n",
" (49.7, 100.0),\n",
" (75.4, 20.0),\n",
" (101.0, 43.74),\n",
" (126.6, 100.0),\n",
" (152.1, 87.89),\n",
" (177.7, 74.21),\n",
" (203.3, 89.23),\n",
" (228.7, 97.07),\n",
" (344.5, 90.7),\n",
" (370.0, 96.37),\n",
" (539.3, 89.12),\n",
" (565.0, 79.88),\n",
" (592.2, 73.09),\n",
" (617.9, 44.03),\n",
" (643.6, 50.58),\n",
" (669.4, 35.55),\n",
" (695.3, 63.8),\n",
" (720.9, 50.0),\n",
" (746.7, 30.0),\n",
" (776.9, 20.0),\n",
" ],\n",
" \"bt\": [\n",
" (0.0, 177.67),\n",
" (15.0, 91.93),\n",
" (30.0, 65.63),\n",
" (45.0, 62.07),\n",
" (60.0, 66.05),\n",
" (75.0, 73.85),\n",
" (90.0, 80.54),\n",
" (105.0, 86.71),\n",
" (120.0, 92.65),\n",
" (135.0, 99.39),\n",
" (150.0, 105.76),\n",
" (165.0, 111.53),\n",
" (180.0, 116.95),\n",
" (195.0, 121.48),\n",
" (210.0, 126.03),\n",
" (225.0, 130.09),\n",
" (240.0, 133.98),\n",
" (255.0, 137.62),\n",
" (270.0, 141.2),\n",
" (285.0, 144.54),\n",
" (300.0, 147.69),\n",
" (315.0, 150.67),\n",
" (330.0, 153.66),\n",
" (345.0, 156.75),\n",
" (360.0, 159.4),\n",
" (375.0, 162.06),\n",
" (390.0, 164.67),\n",
" (405.0, 167.23),\n",
" (420.0, 169.76),\n",
" (435.0, 172.28),\n",
" (450.0, 174.93),\n",
" (465.0, 177.29),\n",
" (480.0, 179.65),\n",
" (495.0, 181.88),\n",
" (510.0, 184.3),\n",
" (525.0, 186.78),\n",
" (540.0, 189.24),\n",
" (555.0, 191.47),\n",
" (570.0, 193.45),\n",
" (585.0, 195.48),\n",
" (600.0, 197.37),\n",
" (615.0, 199.24),\n",
" (630.0, 200.62),\n",
" (645.0, 202.2),\n",
" (660.0, 203.73),\n",
" (675.0, 205.04),\n",
" (690.0, 206.3),\n",
" (705.0, 207.66),\n",
" (720.0, 208.93),\n",
" (735.0, 210.06),\n",
" (750.0, 211.06),\n",
" (765.0, 211.84),\n",
" (780.0, 212.44),\n",
" (795.0, 212.94),\n",
" (810.0, 213.7),\n",
" (825.0, 214.42),\n",
" (840.0, 215.32),\n",
" (851.0, 215.89),\n",
" ],\n",
" },\n",
" # batch 5: 13.98 min, 10 gas steps\n",
" {\n",
" \"mach0\": 233.21,\n",
" \"stop\": 838.71,\n",
" \"turn\": 40.97,\n",
" \"gas\": [\n",
" (0.0, 20.0),\n",
" (60.8, 25.96),\n",
" (86.6, 90.0),\n",
" (564.4, 84.97),\n",
" (620.9, 70.0),\n",
" (684.6, 63.85),\n",
" (710.3, 50.0),\n",
" (768.1, 42.51),\n",
" (794.0, 34.97),\n",
" (820.5, 20.0),\n",
" ],\n",
" \"bt\": [\n",
" (0.0, 179.5),\n",
" (15.0, 91.11),\n",
" (30.0, 65.98),\n",
" (45.0, 62.81),\n",
" (60.0, 66.53),\n",
" (75.0, 72.43),\n",
" (90.0, 80.0),\n",
" (105.0, 87.43),\n",
" (120.0, 94.44),\n",
" (135.0, 101.1),\n",
" (150.0, 106.88),\n",
" (165.0, 112.48),\n",
" (180.0, 117.62),\n",
" (195.0, 122.45),\n",
" (210.0, 126.76),\n",
" (225.0, 130.85),\n",
" (240.0, 134.56),\n",
" (255.0, 138.03),\n",
" (270.0, 141.31),\n",
" (285.0, 144.22),\n",
" (300.0, 147.36),\n",
" (315.0, 150.19),\n",
" (330.0, 153.1),\n",
" (345.0, 156.08),\n",
" (360.0, 158.65),\n",
" (375.0, 161.15),\n",
" (390.0, 163.95),\n",
" (405.0, 166.43),\n",
" (420.0, 168.89),\n",
" (435.0, 171.46),\n",
" (450.0, 173.67),\n",
" (465.0, 175.98),\n",
" (480.0, 178.22),\n",
" (495.0, 180.24),\n",
" (510.0, 182.14),\n",
" (525.0, 184.6),\n",
" (540.0, 186.66),\n",
" (555.0, 188.87),\n",
" (570.0, 191.04),\n",
" (585.0, 193.11),\n",
" (600.0, 195.21),\n",
" (615.0, 197.37),\n",
" (630.0, 199.72),\n",
" (645.0, 201.54),\n",
" (660.0, 203.38),\n",
" (675.0, 205.2),\n",
" (690.0, 206.67),\n",
" (705.0, 207.69),\n",
" (720.0, 208.52),\n",
" (735.0, 209.24),\n",
" (750.0, 209.85),\n",
" (765.0, 210.72),\n",
" (780.0, 211.59),\n",
" (795.0, 212.48),\n",
" (810.0, 213.45),\n",
" (825.0, 214.28),\n",
" (838.7, 215.11),\n",
" ],\n",
" },\n",
"]\n",
"\n",
"# total points: 290\n",
"\n",
"\n",
"print(f\"{len(BATCHES)} batches, {sum(len(b['bt']) for b in BATCHES)} points\")\n",
"print(\"each carries its own gas trace, its own probe reading at charge, and its own\")\n",
"print(\"machine temperature at charge -- the three things that differ between them\")\n",
"for n, b in enumerate(BATCHES, 1):\n",
" print(\n",
" f\" batch {n}: {b['stop'] / 60:5.2f} min, {len(b['gas'])} gas steps, \"\n",
" f\"turn {b['turn']:4.1f} s, drop {b['bt'][-1][1]:6.2f} degC\"\n",
" )"
]
},
{
"cell_type": "markdown",
"id": "a72e11c4",
"metadata": {},
"source": [
"### The turning point is now something the model produces\n",
"\n",
"This is what the third state buys.\n",
"\n",
"At charge the probe is not at bean temperature and not at machine temperature. It has been\n",
"hanging in the recirculating air at 217 °C — and on a Loring that really is what it reads,\n",
"moving air rather than a hot wall — and then fifteen kilos of room-temperature coffee bury\n",
"it. Over the next forty seconds it dives to about 62 °C and turns back up.\n",
"\n",
"Roasters call that minimum the **turning point** and read a lot into it. It is a property of\n",
"the thermocouple as much as of the coffee: the probe is chasing a bean temperature that is\n",
"itself climbing, and the minimum is the instant the two curves cross.\n",
"\n",
"Previously this notebook drew a grey box over the first ninety seconds and said *not\n",
"compared*. With `T_probe` in the model there is nothing to exclude — the dive, the minimum\n",
"and its timing all come out of the integration."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "ec0ef6b8",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"whole curve, 59 samples from the charge -- this page carries the log at 15 s;\n",
"the spec checks the same model against all 558 raw samples, at 1.71 K rms\n",
" rms 1.67 K, worst 7.16 K\n",
"after 120 s, 51 samples:\n",
" rms 1.44 K, worst 2.85 K\n",
"turn: model 63.4 degC at 47 s, measured 64.5 degC at 45 s\n"
]
}
],
"source": [
"b1 = BATCHES[0]\n",
"m_t = np.array([p[0] for p in b1[\"bt\"]])\n",
"m_bt = np.array([p[1] for p in b1[\"bt\"]])\n",
"\n",
"model = np.array([result.at(float(x))[\"probe\"] for x in m_t])\n",
"residual = model - m_bt\n",
"\n",
"fig, (top, bottom) = plt.subplots(\n",
" 2, 1, figsize=(8.4, 4.8), sharex=True, height_ratios=[2.4, 1], gridspec_kw={\"hspace\": 0.15}\n",
")\n",
"\n",
"top.plot(m_t, m_bt, color=ORANGE)\n",
"top.plot(t, probe, color=BLUE)\n",
"label(top, m_t[-1], m_bt[-1], \"measured\", ORANGE, dy=9)\n",
"label(top, t[-1], probe[-1], \"model\", BLUE, dy=-9)\n",
"\n",
"turn = int(np.argmin(probe))\n",
"top.plot([t[turn]], [probe[turn]], marker=\"o\", markersize=8, color=BLUE, zorder=5)\n",
"top.annotate(\n",
" f\"turning point\\n{probe[turn]:.0f} degC at {t[turn]:.0f} s\",\n",
" xy=(t[turn], probe[turn]),\n",
" xytext=(150, 105),\n",
" fontsize=8.5,\n",
" color=MUTED,\n",
" arrowprops={\"arrowstyle\": \"-\", \"color\": MUTED, \"linewidth\": 0.8},\n",
")\n",
"top.set_title(\"Batch 1: the model against the machine, from the charge\")\n",
"top.set_ylabel(\"degC\")\n",
"top.set_xlim(0, STOP * 1.09)\n",
"\n",
"bottom.axhline(0, color=LINE, linewidth=0.8)\n",
"bottom.plot(m_t, residual, color=BLUE)\n",
"bottom.set_ylabel(\"model - measured (K)\")\n",
"minutes(bottom)\n",
"\n",
"fig.legend(\n",
" handles=[\n",
" mpl.lines.Line2D([], [], color=ORANGE, label=\"measured (bean probe)\"),\n",
" mpl.lines.Line2D([], [], color=BLUE, label=\"model (T_probe)\"),\n",
" ],\n",
" loc=\"lower center\",\n",
" ncol=2,\n",
" bbox_to_anchor=(0.5, -0.06),\n",
")\n",
"plt.show()\n",
"\n",
"late = m_t >= SETTLED\n",
"print(f\"whole curve, {len(m_t)} samples from the charge -- this page carries the log at 15 s;\")\n",
"print(\"the spec checks the same model against all 558 raw samples, at 1.67 K rms\")\n",
"print(f\" rms {np.sqrt((residual**2).mean()):.2f} K, worst {np.abs(residual).max():.2f} K\")\n",
"print(f\"after {SETTLED:.0f} s, {late.sum()} samples:\")\n",
"print(\n",
" f\" rms {np.sqrt((residual[late] ** 2).mean()):.2f} K, \"\n",
" f\"worst {np.abs(residual[late]).max():.2f} K\"\n",
")\n",
"print(\n",
" f\"turn: model {probe[turn]:.1f} degC at {t[turn]:.0f} s, \"\n",
" f\"measured {m_bt.min():.1f} degC at {m_t[int(np.argmin(m_bt))]:.0f} s\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "11c0aa30",
"metadata": {},
"source": [
"## What the roaster is reading is late\n",
"\n",
"Nobody roasting coffee watches the temperature. They watch its slope — *rate of rise*, in\n",
"degrees per minute — because that is the thing that has to come down smoothly.\n",
"\n",
"Now that the probe is in the model, there are two rates of rise to compare: the coffee's,\n",
"and the one the roaster can actually see. They are not the same, and the gap between the\n",
"temperatures behind them is not small. Around two minutes the coffee is running **8 K\n",
"ahead** of the number on the display, closing to under a degree by the drop."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "4d223118",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" t = 2.00 min coffee 103.7 probe 95.5 lag +8.22 K\n",
" t = 5.00 min coffee 152.5 probe 149.1 lag +3.49 K\n",
" t = 10.00 min coffee 195.4 probe 193.3 lag +2.11 K\n",
" t = 14.39 min coffee 213.7 probe 213.5 lag +0.16 K\n"
]
}
],
"source": [
"def rate_of_rise(times, values, after=0.0):\n",
" \"\"\"Degrees per minute, and the instants it is defined at.\n",
"\n",
" An event puts two samples at the same instant -- the value on each side of the\n",
" discontinuity, which is the whole point of restarting there -- and a centred difference\n",
" over a zero-width interval is a division by zero. So the duplicates come out first.\n",
" \"\"\"\n",
" keep = np.concatenate(([True], np.diff(times) > 0)) & (times >= after)\n",
" kept, kept_values = times[keep], values[keep]\n",
" return kept, np.gradient(kept_values, kept) * 60.0\n",
"\n",
"\n",
"lag = result.values[result.names.index(\"lag\")]\n",
"\n",
"fig, (top, bottom) = plt.subplots(\n",
" 2, 1, figsize=(8.4, 4.8), sharex=True, height_ratios=[2, 1], gridspec_kw={\"hspace\": 0.15}\n",
")\n",
"\n",
"m_ror_t, m_ror = rate_of_rise(m_t, m_bt, after=SETTLED)\n",
"p_ror_t, p_ror = rate_of_rise(t, probe, after=SETTLED)\n",
"b_ror_t, b_ror = rate_of_rise(t, bean, after=SETTLED)\n",
"top.plot(m_ror_t, m_ror, color=ORANGE, linewidth=1.4, alpha=0.9)\n",
"top.plot(p_ror_t, p_ror, color=BLUE)\n",
"top.plot(b_ror_t, b_ror, color=AQUA, linewidth=1.5, linestyle=(0, (5, 2)))\n",
"label(top, m_ror_t[-1], m_ror[-1], \"measured\", ORANGE)\n",
"top.set_title(\"Rate of rise: what is happening, and what the display shows\")\n",
"top.set_ylabel(\"degC / min\")\n",
"top.set_ylim(0, None)\n",
"top.set_xlim(0, STOP * 1.09)\n",
"\n",
"keep = t >= SETTLED\n",
"bottom.axhline(0, color=LINE, linewidth=0.8)\n",
"bottom.fill_between(t[keep], 0, lag[keep], color=BLUE, alpha=0.16, linewidth=0)\n",
"bottom.plot(t[keep], lag[keep], color=BLUE)\n",
"bottom.set_ylabel(\"coffee - probe (K)\")\n",
"minutes(bottom)\n",
"\n",
"fig.legend(\n",
" handles=[\n",
" mpl.lines.Line2D([], [], color=ORANGE, label=\"measured\"),\n",
" mpl.lines.Line2D([], [], color=BLUE, label=\"model, at the probe\"),\n",
" mpl.lines.Line2D([], [], color=AQUA, linestyle=(0, (5, 2)), label=\"model, the coffee\"),\n",
" ],\n",
" loc=\"lower center\",\n",
" ncol=3,\n",
" bbox_to_anchor=(0.5, -0.06),\n",
")\n",
"plt.show()\n",
"\n",
"for x in [120.0, 300.0, 600.0, STOP]:\n",
" s = result.at(x)\n",
" print(\n",
" f\" t = {x / 60:5.2f} min coffee {s['bean']:6.1f} probe {s['probe']:6.1f} \"\n",
" f\"lag {s['lag']:+5.2f} K\"\n",
" )"
]
},
{
"cell_type": "markdown",
"id": "47a15a51",
"metadata": {},
"source": [
"\n",
"## One parameter set, five batches\n",
"\n",
"Seven numbers describe the machine, and they were fitted against **all five batches at\n",
"once** — 2791 samples — rather than against batch 1 alone. Each batch then keeps three\n",
"things of its own: its recorded gas trace, its probe reading at the charge mark, and its\n",
"machine temperature at the charge mark.\n",
"\n",
"That third one matters more than it looks, and this notebook got it wrong until recently.\n",
"The prose claimed it while `rebound()` rebound only the probe, so every batch was silently\n",
"wearing batch 1's preheat — and the residual grew monotonically down the session, about\n",
"1.5 K per batch, which is the signature of a state that differs between runs and is being\n",
"forced to be the same. `rusty-replicant-32c43` spotted the contradiction between the\n",
"sentence and the code.\n",
"\n",
"Batch 1 followed a 35-minute cold-start preheat; batches 2–5 came off a purge-and-reheat\n",
"between roasts, where the Loring opens a gate, pulls room air until the drum is down around\n",
"77 °C, then closes it and comes back to charge temperature. So the machine really does start\n",
"each batch somewhere slightly different, and it is not a drift — it is one cold start\n",
"followed by four nearly identical turnarounds.\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "5d30adac",
"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 rebound(model: FlatModel, steps, probe0: float, mach0: float) -> FlatModel:\n",
" \"\"\"The same model on a different gas trace, from a different pair of initial readings.\n",
"\n",
" The schedule lives in the equation for `burner` as a chain of comparisons against\n",
" `time`, so swapping it means rebuilding that chain -- which is what the Python\n",
" front-end does when it folds a table, and what `unparse` prints back as\n",
" `if ... elseif ... else`.\n",
"\n",
" Level `i` holds from `steps[i][0]` until `steps[i + 1][0]`, so each comparison carries\n",
" the level *before* it, ending on the last level, which has nothing after it to stop it.\n",
" \"\"\"\n",
" chain: ir.Expr = ir.RealLiteral(float(steps[-1][1]))\n",
" for index in range(len(steps) - 1, 0, -1):\n",
" chain = ir.IfExpr(\n",
" ir.Binary(\"<\", ir.Ref((\"time\",)), ir.RealLiteral(float(steps[index][0]))),\n",
" ir.RealLiteral(float(steps[index - 1][1])),\n",
" chain,\n",
" )\n",
" starts = {\"T_probe\": probe0 + 273.15, \"T_mach\": mach0 + 273.15}\n",
" return replace(\n",
" model,\n",
" equations=tuple(\n",
" replace(eq, rhs=chain)\n",
" if isinstance(eq.lhs, ir.Ref) and eq.lhs.parts == (\"burner\",)\n",
" else eq\n",
" for eq in model.equations\n",
" ),\n",
" variables=tuple(\n",
" replace(v, start=ir.RealLiteral(starts[v.name])) if v.name in starts else v\n",
" for v in model.variables\n",
" ),\n",
" )\n",
"\n",
"\n",
"def roast(batch):\n",
" \"\"\"Run the model against one recorded batch: its gas, its probe, its machine.\"\"\"\n",
" return simulate(\n",
" causalize(rebound(flat, batch[\"gas\"], batch[\"bt\"][0][1], batch[\"mach0\"])),\n",
" stop=batch[\"stop\"],\n",
" tolerance=1e-10,\n",
" points=700,\n",
" )\n",
"\n",
"\n",
"runs = [roast(b) for b in BATCHES]"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "f01743c6",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" batch rms worst after 2 min drop meas model\n",
" 1 1.65K 7.16K 1.42K 216.1 213.6\n",
" 2 2.47K 3.90K 2.62K 216.9 218.7\n",
" 3 2.60K 4.53K 2.75K 215.7 217.2\n",
" 4 3.10K 4.23K 3.29K 215.9 217.7\n",
" 5 3.31K 5.53K 3.52K 215.1 220.1\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Note what this page cannot show you. It carries the log downsampled to 15 s, and the\n",
"worst disagreement on batches 2-4 is in the first second and a half -- where the probe\n",
"is still swinging from the previous drop when the charge mark lands, so its start value\n",
"is contaminated. Against the full 1.5 s log those batches peak at 9 to 11 K there, and\n",
"under 7 K everywhere after. Batch 1 is the only one preceded by a preheat rather than a\n",
"drop, which is why the example carries it.\n"
]
}
],
"source": [
"fig, axes = plt.subplots(1, 5, figsize=(9.6, 2.9), sharex=True, sharey=True)\n",
"\n",
"print(\n",
" f\"{'batch':>6} {'T_mach(0)':>10} {'rms':>8} {'after 2 min':>12} {'drop meas':>10} {'model':>8}\"\n",
")\n",
"for n, (ax, b, run) in enumerate(zip(axes, BATCHES, runs, strict=True), 1):\n",
" bt_t = np.array([p[0] for p in b[\"bt\"]])\n",
" bt_v = np.array([p[1] for p in b[\"bt\"]])\n",
" sim = np.array([run.at(float(x))[\"probe\"] for x in bt_t])\n",
" err = sim - bt_v\n",
" late = bt_t >= SETTLED\n",
"\n",
" ax.plot(bt_t, bt_v, color=ORANGE, linewidth=1.6)\n",
" ax.plot(run.time, run.values[run.names.index(\"probe\")], color=BLUE, linewidth=1.6)\n",
" ax.set_title(f\"batch {n}\", fontsize=10)\n",
" ax.set_xticks([0, 300, 600, 900])\n",
" ax.xaxis.set_major_formatter(lambda s, _: f\"{int(s // 60)}\")\n",
" ax.set_xlabel(\"min\")\n",
" if n == 1:\n",
" ax.set_ylabel(\"degC\")\n",
"\n",
" print(\n",
" f\"{n:>6} {b['mach0']:9.1f}C {np.sqrt((err**2).mean()):7.2f}K \"\n",
" f\"{np.sqrt((err[late] ** 2).mean()):11.2f}K {bt_v[-1]:10.1f} {sim[-1]:8.1f}\"\n",
" )\n",
"\n",
"axes[0].legend(\n",
" handles=[\n",
" mpl.lines.Line2D([], [], color=ORANGE, label=\"measured\"),\n",
" mpl.lines.Line2D([], [], color=BLUE, label=\"model\"),\n",
" ],\n",
" loc=\"lower right\",\n",
" fontsize=8,\n",
")\n",
"fig.tight_layout()\n",
"fig.subplots_adjust(top=0.78)\n",
"fig.suptitle(\n",
" \"Fitted on all five at once, each on its own gas\", x=0.01, y=0.98, ha=\"left\", fontsize=11\n",
")\n",
"plt.show()\n",
"\n",
"print()\n",
"print(\"The residual no longer grows down the session -- and batch 1 is now the *worst* of\")\n",
"print(\"the five rather than the best, which is what you would expect once every batch stops\")\n",
"print(\"wearing its preheat. This page carries the log at 15 s, so it cannot show the one\")\n",
"print(\"place the model still misses badly: the first second or two of batches 2-4, where\")\n",
"print(\"the probe is still swinging from the previous drop when the charge mark lands.\")"
]
},
{
"cell_type": "markdown",
"id": "bb441586",
"metadata": {},
"source": [
"\n",
"## The knobs a roaster actually has\n",
"\n",
"Four of them, and all four are parameters of the flat model, so tuning them is\n",
"`dataclasses.replace` and a re-causalization — the same trick as\n",
"[Turning the knobs](exploring.ipynb).\n",
"\n",
"* **charge** is `m_green`. The measured charge is 13.8 kg, not the machine's rated 15, and\n",
" because `m_dry` and the starting water are computed from it, changing it moves the heat\n",
" capacity *and* the mass balance.\n",
"* **moisture** is the fraction of that charge which is water.\n",
"* **preheat** is where the machine starts — the number that separates batch 1 from the\n",
" four that follow it.\n",
"* **gas trim** scales `Q_max`, standing in for a different line pressure or a dirtier\n",
" burner.\n",
"\n",
"The readout is what a roaster cares about: how hot it was when it came out, how long it took\n",
"to reach 196 °C — about where first crack starts, and the recorded batch hit it at 11:25 —\n",
"and how much of the charge left as steam.\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "82e1cd65",
"metadata": {},
"outputs": [],
"source": [
"def tuned(model: FlatModel, values: dict, preheat: float | None = None) -> FlatModel:\n",
" \"\"\"A copy with those parameters rebound, and optionally a different preheat.\n",
"\n",
" `m_dry` and the initial water are *bindings over other parameters* rather than\n",
" literals, so rebinding `m_green` or `moisture` carries through to both of them without\n",
" anything here having to know that it should.\n",
" \"\"\"\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(\n",
" v,\n",
" binding=ir.RealLiteral(float(values[v.name])) if v.name in values else v.binding,\n",
" start=(\n",
" ir.RealLiteral(preheat + 273.15)\n",
" if preheat is not None and v.name == \"T_mach\"\n",
" else v.start\n",
" ),\n",
" )\n",
" for v in model.variables\n",
" ),\n",
" )\n",
"\n",
"\n",
"def crack(run, threshold: float = 196.0) -> float | None:\n",
" \"\"\"When the probe first reads `threshold`, interpolated between output points.\"\"\"\n",
" times = run.time\n",
" reading = run.values[run.names.index(\"probe\")]\n",
" hit = np.flatnonzero(reading >= threshold)\n",
" if not len(hit):\n",
" return None\n",
" k = int(hit[0])\n",
" if k == 0:\n",
" return float(times[0])\n",
" lo, hi = reading[k - 1], reading[k]\n",
" return float(times[k - 1] + (threshold - lo) / (hi - lo) * (times[k] - times[k - 1]))"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "b8d3257c",
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "be16838792a64de2a4f6a36e825bbdb9",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"interactive(children=(FloatSlider(value=15.0, description='charge (kg)', max=18.0, min=8.0, step=0.5), FloatSl…"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
" None>"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import ipywidgets as widgets\n",
"\n",
"\n",
"def show(charge: float, moisture: float, preheat: float, trim: float) -> None:\n",
" run = simulate(\n",
" causalize(\n",
" tuned(\n",
" flat,\n",
" {\"m_green\": charge, \"moisture\": moisture, \"Q_max\": 1.32e5 * trim},\n",
" preheat=preheat,\n",
" )\n",
" ),\n",
" stop=STOP,\n",
" tolerance=1e-9,\n",
" points=700,\n",
" )\n",
" reading = run.values[run.names.index(\"probe\")]\n",
" fc = crack(run)\n",
"\n",
" fig, ax = plt.subplots(figsize=(8.4, 3.4))\n",
" ax.plot(m_t, m_bt, color=ORANGE, linewidth=1.4)\n",
" ax.plot(run.time, reading, color=BLUE)\n",
" label(ax, m_t[-1], m_bt[-1], \"as roasted\", ORANGE, dy=9)\n",
" label(ax, run.time[-1], reading[-1], \"this batch\", BLUE, dy=-9)\n",
" if fc is not None:\n",
" ax.plot([fc], [196.0], marker=\"o\", markersize=8, color=BLUE, zorder=5)\n",
" ax.axhline(196.0, color=LINE, linewidth=0.8, zorder=0)\n",
"\n",
" end = run.at(STOP)\n",
" when = f\"{fc / 60:.2f} min\" if fc is not None else \"never\"\n",
" ax.set_title(\n",
" f\"drop {end['probe']:.1f} degC · 196 degC at {when} · \"\n",
" f\"{end['loss']:.2f} % lost as steam\"\n",
" )\n",
" ax.set_ylabel(\"degC\")\n",
" ax.set_ylim(0, 260)\n",
" ax.set_xlim(0, STOP * 1.09)\n",
" minutes(ax)\n",
" ax.legend(\n",
" handles=[\n",
" mpl.lines.Line2D([], [], color=ORANGE, label=\"batch 1, as roasted\"),\n",
" mpl.lines.Line2D([], [], color=BLUE, label=\"model, with these knobs\"),\n",
" ],\n",
" loc=\"lower right\",\n",
" )\n",
" plt.show()\n",
"\n",
"\n",
"widgets.interact(\n",
" show,\n",
" charge=widgets.FloatSlider(value=13.8, min=8.0, max=18.0, step=0.2, description=\"charge (kg)\"),\n",
" moisture=widgets.FloatSlider(\n",
" value=0.13, min=0.07, max=0.18, step=0.005, readout_format=\".3f\", description=\"moisture\"\n",
" ),\n",
" preheat=widgets.FloatSlider(\n",
" value=237.4, min=200.0, max=270.0, step=1.0, description=\"preheat (degC)\"\n",
" ),\n",
" trim=widgets.FloatSlider(value=1.0, min=0.7, max=1.3, step=0.02, description=\"gas trim\"),\n",
")"
]
},
{
"cell_type": "markdown",
"id": "9cbb77ba",
"metadata": {},
"source": [
"\n",
"## What the weighing changed\n",
"\n",
"When this notebook was first written it ended with a paragraph asking for one number: the\n",
"weight of a batch, green in and roasted out. It arrived — off a screen recording of the\n",
"roaster's own session, same four batches:\n",
"\n",
"| | green | roasted | lost |\n",
"|---|---|---|---|\n",
"| PR-1913 | 13.8 kg | 11.76 kg | 14.78 % |\n",
"| PR-1914 | 13.8 kg | 11.74 kg | 14.93 % |\n",
"| PR-1915 | 13.8 kg | 11.76 kg | 14.78 % |\n",
"| PR-1916 | 13.8 kg | 11.73 kg | 15.00 % |\n",
"\n",
"Two things fell out of it, and the second is the interesting one.\n",
"\n",
"**The charge was 13.8 kg, not 15.** The model had been using the machine's *rated* capacity\n",
"because the bridge capture carried no weight — a rated capacity is not a charge, and that\n",
"alone was an 8 % error in the heat capacity of everything.\n",
"\n",
"**The drying law was wrong, not merely mistuned.** The old model gated evaporation at\n",
"100 °C, on the reasoning that water leaves when water boils, and removed 2.7 % of the charge\n",
"where ~11 points of water should go. Forced to remove the right mass it fell apart entirely.\n",
"Sweeping the onset downward fixes the fit and the mass balance *together*, monotonically,\n",
"all the way down to ambient:\n",
"\n",
"| onset | curve rms | water lost | dried by half-time |\n",
"|---|---|---|---|\n",
"| 100 °C | 1.806 K | 11.24 % | 38.7 % |\n",
"| 60 °C | 1.492 K | 11.27 % | 51.9 % |\n",
"| 40 °C | 1.373 K | 11.37 % | 58.1 % |\n",
"| 20 °C | 1.232 K | 11.41 % | 64.3 % |\n",
"\n",
"A coffee bean is not a puddle. It is porous, and it gives up bound water by diffusion the\n",
"whole time it is warm — which is why roasters call the *first* third of a roast the drying\n",
"phase. `T_dry = 40 °C` is a modelling choice costing 0.14 K against no threshold at all; it\n",
"keeps the crossing a real event instead of one that fires at `t = 0` and means nothing.\n",
"\n",
"**How this was nearly missed.** Sweeping the drying with the seven machine parameters held\n",
"fixed, *no* combination reproduces both the curve and the mass — and the obvious reading is\n",
"that a term is missing, with first crack being genuinely exothermic making that a tempting\n",
"one to reach for. But those seven had been fitted against a model whose latent sink was four\n",
"times too small, so `UA`, `G_loss` and `Q_max` had already absorbed the missing heat. Adding\n",
"it back on top of them cannot work. Refitting everything together is what showed the shape\n",
"had been fine all along.\n",
"\n",
"## What is still assumed\n",
"\n",
"`moisture` is still a parameter, not a measurement. The weighing gives *total* loss, and a\n",
"roast also sheds dry matter — CO2, volatiles, chaff, three or four points of it — which this\n",
"model has no term for. So the water target is the measured 14.87 % less an assumed 3.5 %\n",
"organic loss, and that 3.5 % is now the softest number in the whole example. A green\n",
"moisture meter would turn both it and `moisture` into one measurement and one residual.\n",
"\n",
"The rest, and the measured points these figures are drawn against, is in\n",
"[`RoastingMachine.toml`](https://gitlab.com/jorgeecardona/pymodelica/-/blob/main/examples/models/RoastingMachine.toml).\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"pygments_lexer": "ipython3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}