01 - A house with a heat pump: the reference implementation¶
Learning objectives
- Follow the pipeline from the Trano YAML to the CasADi model: what the FastModel of a house contains, which inputs the controller decides and which ones it has to forecast, where the energy is metered.
- Read an economic MPC plan: how the comfort band, the price signal and the horizon shape the heating of a heat pump house with a buffer tank and a radiator.
- Run the closed loop against a building that is not the control model (perturbed envelope, noisy sensors) and compare it with the controllers written in the YAML.
- See why state estimation matters: the controller needs the walls, the water and the buffer, it measures the air temperature only.
Everything here calls the library; the closed loops come from tutorials/cache/house/
(uv run python tutorials/scripts/build_cache.py house). The plan widget solves live.
Runtimes (hourly steps, this container): Trano and rumoca about 10 s per model; one 24-hour plan 5 to 20 s; one two-day closed loop with moving horizon estimation about ten minutes; the rule-based reference a second.
import sys
if "google.colab" in sys.modules: # Colab: fetch the tutorials, install mpctrano and rumoca (a no-op elsewhere)
import subprocess
subprocess.run( # a fixed command; the clone is skipped when it exists
["git", "clone", "-q", "--depth", "1", "https://github.com/andoludo/mpctrano", "/content/mpctrano"], check=False
)
sys.path.insert(0, "/content/mpctrano")
from tutorials import _colab
_colab.setup()
import pandas as pd
from IPython.display import display
from mpctrano import ComfortSpec, ModelPredictiveController, MPCSettings
from mpctrano.notebook import closed_loop_explorer, kpi_table, plan_explorer, quick_mode
from mpctrano.plot import plot_closed_loop, plot_plan, plot_states
from tutorials import _house
pd.set_option("display.width", 200)
print(f"quick mode: {quick_mode()}; cached runs: {_house.DAYS:g} days from {_house.START.date()}")
quick mode: False; cached runs: 2 days from 2023-01-16
1. From the YAML to the CasADi model¶
tutorials/buildings/house/base.yaml describes a 100 m2 house: cavity walls, triple glazing, a concrete slab,
a 3.5 kW air-to-water heat pump with a 150 l buffer tank, a mixing valve and a low-temperature radiator behind
a thermostatic valve. Trano renders it with its FastModelMPC library (physics close to the Buildings library:
layered walls, window optics, water nodes) where the controls of the systems are inputs, and with its
FastModel library where the controllers of the YAML close the loop (the rule-based reference). rumoca
translates each into an explicit CasADi ODE xdot = f(t, x, u, d, p).
house = _house.building("base")
model = house.model
print(f"{model.n_states} states, {model.n_controls} controls, {model.n_disturbances} disturbances, "
f"{model.n_parameters} parameters") # fmt: skip
print("controls ", model.control_names)
print("disturbances", model.disturbance_names)
print("comfort ", model.comfort_states)
print("meters ", model.meters)
print("states ", model.state_names)
27 states, 3 controls, 9 disturbances, 581 parameters
controls ['boiler_001_y.y', 'three_way_valve_001_y.y', 'valve_001_y.y']
disturbances ['weather_inputs.TDryBul', 'weather_inputs.TBlaSky', 'weather_inputs.winSpe', 'weather_inputs.winDir', 'weather_inputs.HDirNor', 'weather_inputs.HDifHor', 'weather_inputs.HGloHor', 'occupancy_1_gai.QCon', 'occupancy_1_gai.QRad']
comfort ['space_001.TAir']
meters {'boiler_001.EEle': 'electricity'}
states ['space_001.TAir', 'externalwall_0.T[1]', 'externalwall_0.T[2]', 'externalwall_0.T[3]', 'externalwall_0.T[4]', 'externalwall_0.T[5]', 'externalwall_1.T[1]', 'externalwall_1.T[2]', 'externalwall_1.T[3]', 'externalwall_1.T[4]', 'externalwall_1.T[5]', 'externalwall_2.T[1]', 'externalwall_2.T[2]', 'externalwall_2.T[3]', 'externalwall_2.T[4]', 'externalwall_2.T[5]', 'window_0.T[1]', 'window_0.T[2]', 'window_0.T[3]', 'window_1.T[1]', 'window_1.T[2]', 'window_1.T[3]', 'flooronground_0.T[1]', 'flooronground_0.T[2]', 'radiator_001.T[1]', 'boiler_001.T', 'boiler_001.EEle']
The roles come from the Modelica model itself: an input with a minimum and a maximum is a control (the
compressor speed boiler_001_y, the share of supply water of the mixing valve, the opening of the radiator
valve, all between 0 and 1), any other input a disturbance (seven weather channels and the internal gains
of the occupancy), the air temperature TAir is the comfort variable and the state in kWh is the electricity
meter of the heat pump. The controller never sees a COP: it optimises the valves and the compressor, and the
electricity follows from the heat pump model and its meter.
A FastModel starts with every temperature at its start value (walls at the set point, buffer at the supply
temperature). Every run of the tutorials starts from the state the controllers of the YAML reach after three
days (spin_up): the night setback has cooled the house at the start of the period.
x0 = pd.Series(house.initial_state(), index=model.state_names)
display((x0[[name for name in model.state_names if name.endswith(("TAir", ".T", "T[1]"))]] - 273.15).round(1))
forecast = house.forecaster.realisation(_house.START, 48, 3600.0)
display(forecast.describe().round(1).T[["mean", "min", "max"]])
space_001.TAir 16.4 externalwall_0.T[1] 5.9 externalwall_1.T[1] 5.4 externalwall_2.T[1] 5.1 window_0.T[1] 5.6 window_1.T[1] 5.6 flooronground_0.T[1] 12.2 radiator_001.T[1] 17.4 boiler_001.T 36.7 dtype: float64
| mean | min | max | |
|---|---|---|---|
| weather_inputs.TDryBul | 277.7 | 274.5 | 279.3 |
| weather_inputs.TBlaSky | 269.8 | 265.2 | 273.3 |
| weather_inputs.winSpe | 4.8 | 2.5 | 6.5 |
| weather_inputs.winDir | 4.1 | 3.6 | 5.2 |
| weather_inputs.HDirNor | 47.4 | 0.0 | 289.0 |
| weather_inputs.HDifHor | 20.7 | 0.0 | 102.0 |
| weather_inputs.HGloHor | 32.0 | 0.0 | 175.0 |
| occupancy_1_gai.QCon | 2.2 | 0.0 | 3.6 |
| occupancy_1_gai.QRad | 1.2 | 0.0 | 1.9 |
2. One economic MPC plan¶
The controller minimises, over the horizon, the electricity bought (read from the meter of the heat pump) at the tariff price, plus penalties on the comfort band, subject to the ODE collocated on Radau points and the bounds of the controls; IPOPT solves the nonlinear programme. With a flat price it heats just enough and as late as the slow radiator allows; with a time-varying price it charges the buffer, the water and the walls when electricity is cheap. Try it: the widget re-solves the plan (5 to 20 s; the next solves warm start from the previous plan).
tariffs = _house.tariffs()
display(plan_explorer(house, tariffs))
The widget needs a live kernel; the plan of its default selection (day-ahead prices, the standard band, a 24-hour horizon), for the static rendering of this notebook:
horizon = 6 if quick_mode() else 24
controller = ModelPredictiveController(model, MPCSettings(horizon=horizon), ComfortSpec(), tariffs["day-ahead"])
plan = controller.solve(house.initial_state(), house.forecaster.forecast(_house.START, horizon, 3600.0))
display(plot_plan(plan, model, "base house, day-ahead prices"))
What to look for:
- The morning recovery: the setback ends at 7h; the plan starts the compressor before, so that the air is in the band at 7h, and earlier when the night is cheap.
- The buffer is the fast storage: the heat pump charges it at full speed in the cheap hours and the mixing valve doses the supply of the radiator.
- The spike at 17h empties the compressor in the evening: the house coasts on the buffer and the walls.
- The comfort weight trades a few kelvin-hours for cost when the band is tight.
3. Closed loop against the true building¶
A plan is not a result: the controller applies the first step, the building responds, the sensor reports the air temperature, the state is estimated and the plan is rebuilt. The true building has its envelope conductances and capacities, air capacity and infiltration drawn at random (20 %) and a sensor noise of 0.1 K; the controller keeps the Trano model. The runs are cached; Recompute runs one again (about ten minutes).
runs = {label: parts for label, parts in _house.runs().items() if parts["estimator"] == "mhe"}
dimensions = {
"variant": list(_house.VARIANTS),
"price": [*_house.PRICES, "day-ahead + evening spike"],
"band": list(_house.BANDS),
"mismatch": ["mismatch"],
}
display(closed_loop_explorer(_house.CACHE, runs, _house.compute, dimensions))
spike = _house.cached(
variant="base",
price="day-ahead + evening spike",
band="standard",
mismatch="mismatch",
estimator="mhe",
days=_house.DAYS,
)
display(plot_closed_loop(spike, "base house, day-ahead prices with an evening spike (the widget's static fallback)"))
MPC against the controllers of the YAML¶
The same building, period, tariff and comfort band, controlled by the thermostat, heating curve and buffer
controller written in the YAML (the closed-loop FastModel, simulated in a second): this is the reference an
MPC has to beat.
comparison = {}
for variant in _house.VARIANTS:
for price in _house.PRICES:
parts = {"variant": variant, "price": price, "band": "standard", "mismatch": "mismatch"}
comparison[f"{variant} | {price} | MPC"] = _house.cached(**parts, estimator="mhe", days=_house.DAYS)
comparison[f"{variant} | {price} | rule based"] = _house.reference(**parts)
display(kpi_table(comparison, ("total_cost", "energy_kwh", "comfort_violation_kh", "peak_power_kw", "low_price_share")))
display(plot_closed_loop(comparison["base | day-ahead | MPC"], "base house, day-ahead prices, true building"))
| total_cost | energy_kwh | comfort_violation_kh | peak_power_kw | low_price_share | |
|---|---|---|---|---|---|
| base | flat | MPC | 3.17 | 12.49 | 1.54 | 1.14 | 0.00 |
| base | flat | rule based | 3.34 | 13.16 | 47.63 | 0.50 | 0.00 |
| base | day-ahead | MPC | 3.21 | 12.60 | 1.40 | 0.94 | 0.58 |
| base | day-ahead | rule based | 3.48 | 13.16 | 47.63 | 0.50 | 0.40 |
| old | flat | MPC | 6.44 | 25.35 | 2.59 | 2.04 | 0.00 |
| old | flat | rule based | 6.71 | 26.43 | 57.91 | 0.99 | 0.00 |
| old | day-ahead | MPC | 6.46 | 25.64 | 2.53 | 1.41 | 0.59 |
| old | day-ahead | rule based | 6.87 | 26.43 | 57.91 | 0.99 | 0.49 |
The price of not knowing the building¶
The same price and band with the plant equal to the control model (and its true state) and with the true building: the cost difference is the price of the model error, the comfort violations show where the controller was over-confident.
pairs = {}
for price in _house.PRICES:
for mismatch, estimator in (("perfect", "perfect"), ("mismatch", "mhe")):
pairs[f"{price} | {mismatch}"] = _house.cached(
variant="base", price=price, band="standard", mismatch=mismatch, estimator=estimator, days=_house.DAYS
)
display(kpi_table(pairs))
| total_cost | energy_kwh | sold_kwh | comfort_violation_kh | peak_power_kw | low_price_share | failures | |
|---|---|---|---|---|---|---|---|
| flat | perfect | 3.02 | 11.87 | 0.0 | 0.03 | 1.07 | 0.00 | 0.0 |
| flat | mismatch | 3.17 | 12.49 | 0.0 | 1.54 | 1.14 | 0.00 | 0.0 |
| day-ahead | perfect | 3.06 | 11.98 | 0.0 | 0.02 | 0.88 | 0.58 | 0.0 |
| day-ahead | mismatch | 3.21 | 12.60 | 0.0 | 1.40 | 0.94 | 0.58 | 0.0 |
4. Why state estimation matters¶
The controller needs the full state: the wall layers, the water of the radiator and the buffer. It measures the air temperature only. The moving horizon estimator fits the trajectory of the model to the measurements of the last hours; below, its estimate of the buffer and of the inner layer of a wall against the true building, then three runs on the same true building: the true state handed to the controller, moving horizon estimation, and no estimation (the measured air temperature, the other states held at their expected start values).
estimation = {
estimator: _house.cached(
variant="base", price="day-ahead", band="standard", mismatch="mismatch", estimator=estimator, days=_house.DAYS
)
for estimator in ("perfect", "mhe", "none")
}
display(plot_states(estimation["mhe"], ["boiler_001.T", "externalwall_0.T[5]"], "Moving horizon estimation"))
display(kpi_table(estimation, ("total_cost", "energy_kwh", "comfort_violation_kh", "max_violation", "failures")))
| total_cost | energy_kwh | comfort_violation_kh | max_violation | failures | |
|---|---|---|---|---|---|
| perfect | 3.19 | 12.50 | 2.31 | 0.10 | 0.0 |
| mhe | 3.21 | 12.60 | 1.40 | 0.13 | 0.0 |
| none | 2.65 | 9.91 | 40.92 | 2.32 | 0.0 |
Key takeaways¶
- The MPC beats the controllers of the YAML on both counts: over the same two January days, on the true building, it costs 5 to 8 % less (3.21 against 3.48 EUR under day-ahead prices for the base house, 6.46 against 6.87 EUR for the old one) with 1.4 K.h of comfort violation against 48 K.h for the thermostat, which hovers at the bottom of the band; it uses about 5 % less energy because the house is never overheated.
- The price signal moves the heating, not the bill: under day-ahead prices the MPC buys 58 % of its energy in the cheap third of the hours (40 % for the thermostat, by chance of the night setback); with so flat a winter price profile the saving is a few cents, the flexibility is what the next tutorial sells.
- Model error costs about 5 %: against the true building (envelope 20 % off, noisy sensor) the same controller pays 3.17 instead of 3.02 EUR and lets the air drift 1.5 K.h outside the band; with the model as the plant the violation is nil.
- State estimation is not optional: with the air temperature only (the other states held at their expected values) the controller under-heats, saves 17 % and violates the band by 41 K.h with 2.3 K at worst; the moving horizon estimate (1.4 K.h) even beats the true state handed over (2.3 K.h), because its bias absorbs part of the model error.
Further reading¶
tutorials/references.md, section A: Drgoňa et al. (2020) for MPC in buildings, Blum et al. (2021) for the KPIs used here.docs/mpc.mdfor the optimisation problem anddocs/estimation.mdfor the estimator.- Next:
02_flexibility_use_casesquantifies what this flexibility is worth under real price signals, capacity tariffs and explicit activation requests.