03 - Calibration: the YAML declares the parameters, IPOPT finds them¶

Learning objectives

Declare the parameters to calibrate in the Trano YAML itself ({value, parameter: true, min, max}), keep it the single source of truth, and see which parameters of the CasADi model each field moves.
Calibrate the closed-loop model (the controllers of the YAML inside, no control input) on monitoring data a real building would give: a noisy air temperature and an electricity meter, no valve positions.
Read how IPOPT solves it: the three fields as decision variables, the closed-loop model simulated with CVODES inside a weighted least squares with priors, the gradient from the adjoint sensitivities of the integrator.
Compare the model before and after: the temperatures and the energy on the training days and on two validation days, and the recovered values against the true ones.

The true building, its data and the calibration come from tutorials/cache/calibration/ (make cache); the calibration is solved live when the notebook does not run in quick mode (a few minutes).

Runtimes (this container): one regeneration of the model by Trano and rumoca about 15 s (one per declared parameter, for the sensitivities); the calibration of three days of hourly data, 38 states, two to five minutes.

In [1]:

Copied!





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 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()

In [2]:

Copied!





import pandas as pd
from IPython.display import display

from mpctrano import CalibrationSettings, Calibrator
from mpctrano.calibration import evaluate, simulate_fit
from mpctrano.notebook import quick_mode, timed
from mpctrano.plot import plot_fit
from tutorials import _calibration as study

pd.set_option("display.width", 200)
print(f"quick mode: {quick_mode()}")
import pandas as pd
from IPython.display import display

from mpctrano import CalibrationSettings, Calibrator
from mpctrano.calibration import evaluate, simulate_fit
from mpctrano.notebook import quick_mode, timed
from mpctrano.plot import plot_fit
from tutorials import _calibration as study

pd.set_option("display.width", 200)
print(f"quick mode: {quick_mode()}")

quick mode: False

1. The parameters, declared in the YAML¶

tests/models/single_zone_calibration.yaml is the heat pump house of tutorial 01 with its cavity wall spelled out and three fields annotated: the conductivity of the insulation, its thickness and the infiltration rate. CalibrationSpec reads them by path (the convention of trano.runtime) and writes the flat copy Trano reads.

In [3]:

Copied!





spec = study.spec()
text = study.YAML.read_text()
print(text[: text.index("    id: SPACE:001")])
display(pd.DataFrame([parameter.model_dump(by_alias=True) for parameter in spec.parameters]).set_index("path"))
spec = study.spec()
text = study.YAML.read_text()
print(text[: text.index("    id: SPACE:001")])
display(pd.DataFrame([parameter.model_dump(by_alias=True) for parameter in spec.parameters]).set_index("path"))

# The single-zone heat pump house with calibration parameters declared in the YAML: a numeric field written as
# {value, parameter: true, min, max} is calibrated within its bounds (mpctrano.yaml_parameters); Trano reads the
# flat copy mpctrano writes. The cavity wall is spelled out so that its insulation can be calibrated.
default: !include_default
material:
  - id: ROCKWOOL:CAL
    density: 60.0
    longwave_emissivity: 0.9
    shortwave_emissivity: 0.6
    specific_heat_capacity: 840.0
    thermal_conductivity: {value: 0.04, parameter: true, min: 0.025, max: 0.08}
constructions:
  - id: CAVITYWALL:CAL
    layers:
      - material: BRICK:001
        thickness: 0.08
      - material: ROCKWOOL:CAL
        thickness: {value: 0.1, parameter: true, min: 0.04, max: 0.2}
      - material: BRICKHOLLOW:001
        thickness: 0.14
      - material: GYPSUM:001
        thickness: 0.015
spaces:
  - occupancy:
      parameters:
        occupancy: 3600*{7, 22}
        gain: "[40; 75; 40]"
        heat_gain_if_occupied: 1/7/3
    parameters:
      floor_area: 100.0
      average_room_height: 2.5
      ach: {value: 0.3, parameter: true, min: 0.05, max: 1.5}

	value	min	max
path
material[ROCKWOOL:CAL].thermal_conductivity	0.04	0.025	0.08
constructions[CAVITYWALL:CAL].layers[1].thickness	0.10	0.040	0.20
spaces[SPACE:001].parameters.ach	0.30	0.050	1.50

Which parameters of the CasADi model a field moves is measured: the model is regenerated with the field raised by 20 % and the exponent of every parameter that moved is read. The conductivity moves the conductance of the insulation layer of the three walls (exponent 1), the thickness moves it with -1 and the capacities of the neighbouring nodes a little, the infiltration rate is the ACH of the zone.

The table also says what not to calibrate together: the conductivity and the thickness enter the model through the same conductance (k / x, exponents +1 and -1), the capacities they move a little apart; no data on the air temperature tells them apart. The study keeps the thickness as drawn (a drawing is a measurement) and calibrates the conductivity and the infiltration rate.

In [4]:

Copied!





for group in study.declared_groups():
    moved = {
        name: round(exponent, 3)
        for name, exponent in group.exponents.items()
        if not name.startswith(("externalwall_1", "externalwall_2"))
    }
    print(f"{group.name}: value {group.value:g} in [{group.lower:g}, {group.upper:g}] -> {moved} (and the other walls)")
groups = study.groups()
print("calibrated:", [group.name for group in groups])
for group in study.declared_groups():
    moved = {
        name: round(exponent, 3)
        for name, exponent in group.exponents.items()
        if not name.startswith(("externalwall_1", "externalwall_2"))
    }
    print(f"{group.name}: value {group.value:g} in [{group.lower:g}, {group.upper:g}] -> {moved} (and the other walls)")
groups = study.groups()
print("calibrated:", [group.name for group in groups])

material[ROCKWOOL:CAL].thermal_conductivity: value 0.04 in [0.025, 0.08] -> {'externalwall_0.constructionType.k[2]': 1.0, 'externalwall_0.G[2]': 1.0} (and the other walls)
constructions[CAVITYWALL:CAL].layers[1].thickness: value 0.1 in [0.04, 0.2] -> {'externalwall_0.constructionType.x[2]': 1.0, 'externalwall_0.C[2]': 0.043, 'externalwall_0.C[3]': 0.05, 'externalwall_0.G[2]': -1.0} (and the other walls)
spaces[SPACE:001].parameters.ach: value 0.3 in [0.05, 1.5] -> {'space_001.ACH': 1.0} (and the other walls)
calibrated: ['material[ROCKWOOL:CAL].thermal_conductivity', 'spaces[SPACE:001].parameters.ach']

2. The true building and its monitoring data¶

The true building is the same YAML with other values of the three parameters: a worse insulation, thinner than drawn, and a leakier house. It is run by its own controllers (the closed-loop FastModel: thermostat, heating curve, buffer control) over five cold days of January, after two days of spin-up; the monitoring system records the weather, the air temperature with 0.1 K of noise and the electricity meter with 20 Wh of noise. It does not record the valves or the compressor: the model being calibrated has its controllers inside.

In [5]:

Copied!





model = study.model()
print(
    "states",
    model.n_states,
    "| controls",
    model.control_names,
    "| measured",
    model.comfort_states,
    "| meters",
    list(model.meters),
)
truth = pd.Series(spec.nominal, name="assumed (YAML)").to_frame()
truth["true value"] = pd.Series(study.TRUTH).reindex(truth.index).fillna(truth["assumed (YAML)"])
display(truth)
data, true_states = study.recorded()
training, validation = study.split(data)
print(f"{len(training) - 1} training hours from {training.index[0]}, {len(validation) - 1} validation hours")
display(data[[*model.comfort_states, *model.meters]].describe().round(2).T)
model = study.model()
print(
    "states",
    model.n_states,
    "| controls",
    model.control_names,
    "| measured",
    model.comfort_states,
    "| meters",
    list(model.meters),
)
truth = pd.Series(spec.nominal, name="assumed (YAML)").to_frame()
truth["true value"] = pd.Series(study.TRUTH).reindex(truth.index).fillna(truth["assumed (YAML)"])
display(truth)
data, true_states = study.recorded()
training, validation = study.split(data)
print(f"{len(training) - 1} training hours from {training.index[0]}, {len(validation) - 1} validation hours")
display(data[[*model.comfort_states, *model.meters]].describe().round(2).T)

states 43 | controls [] | measured ['space_001.TAir'] | meters ['boiler_001.EEle']

	assumed (YAML)	true value
material[ROCKWOOL:CAL].thermal_conductivity	0.04	0.06
constructions[CAVITYWALL:CAL].layers[1].thickness	0.10	0.10
spaces[SPACE:001].parameters.ach	0.30	0.60

72 training hours from 2023-01-09 00:00:00, 48 validation hours

	count	mean	std	min	25%	50%	75%	max
space_001.TAir	121.0	290.00	0.66	288.53	289.42	290.13	290.46	291.30
boiler_001.EEle	121.0	36.35	13.63	12.01	24.33	38.13	48.42	57.11

3. Before: the model with the values of the YAML¶

An open-loop simulation of the model from the first measurement, with the assumed parameters: the house loses less heat than the true one, so it stays warmer and draws less electricity.

In [6]:

Copied!





before = simulate_fit(model, training)
metrics_before, _ = evaluate(model, training)
display(plot_fit(before, "before calibration: training days (dotted: measured)"))
display(metrics_before.frame().round(3))
before = simulate_fit(model, training)
metrics_before, _ = evaluate(model, training)
display(plot_fit(before, "before calibration: training days (dotted: measured)"))
display(metrics_before.frame().round(3))

	rmse	mae	max_error	r2
space_001	1.351	1.267	2.012	-3.045

4. How IPOPT calibrates it¶

The calibration is one nonlinear programme over the whole training period, solved by single shooting (CalibrationSettings.method, chosen automatically for a closed-loop model):

decision variables: one per declared parameter, log(field / value), bounded by log(min / value) and log(max / value); nothing else. The parameters of the model are expressions of the fields (p = p0 * (field / value) ** e, the exponents measured above; a conductivity and a thickness both move the conductance of the insulation layer, their factors compound).
the simulation inside the objective: CVODES integrates the closed-loop model hour by hour, with the parameters of the candidate, from the start values of the Modelica model through two days of spin-up before the data (so that the walls, the water and the controller states at the start are those of the candidate, not of the nominal model) and then over the training days with their weather.
the gradient: central differences of the objective with a step of 5 % of a field, 2n + 1 simulations of a second each; IPOPT's Hessian is a quasi-Newton (L-BFGS) one and it takes its full steps. A thermostat that flips one integration step earlier makes the objective jumpy at the scale of a tenth of a percent of a field: the adjoint sensitivities of the integrator and small finite differences read that noise and stall IPOPT short of the optimum, a wide step reads the slope of the valley through it and a sufficient-decrease test would refuse steps for the same jumps.
objective: the squared errors of the simulated air temperature against the measurements (weighted by the sensor noise), of the energy metered over each hour against the meter (weighted by the meter noise), and the squared deviations of the fields from their prior (the YAML).

The models with control inputs (the MPC models of the other tutorials) are calibrated on their Radau collocation instead, the states as decision variables and the dynamics as constraints: well conditioned over long periods and poor initial parameters. Not here: the thermostat, the hysteresis of the heat pump controller and the valve limits make the collocated dynamics of a closed-loop model non-smooth, and IPOPT crawls for hundreds of iterations without converging; a simulation follows these switches in under a second, and three parameters need a handful of iterations. The log below is IPOPT's own (one line per iteration: the objective, the dual infeasibility, the step).

In [7]:

Copied!





if quick_mode():
    calibrated, result = study.calibrate()
    print("cached result:", result.status, "in", result.iterations, "iterations,", f"{result.solve_time:.0f} s")
else:
    with timed("calibration"):
        settings = CalibrationSettings(groups=groups, verbose=True, max_iterations=60)
        calibrated, result = Calibrator(model, settings).fit(training, validation)
print(
    f"{result.method}: {result.variables} decision variables, {result.iterations} iterations, "
    f"{result.solve_time:.0f} s, {result.status}"
)
if quick_mode():
    calibrated, result = study.calibrate()
    print("cached result:", result.status, "in", result.iterations, "iterations,", f"{result.solve_time:.0f} s")
else:
    with timed("calibration"):
        settings = CalibrationSettings(groups=groups, verbose=True, max_iterations=60)
        calibrated, result = Calibrator(model, settings).fit(training, validation)
print(
    f"{result.method}: {result.variables} decision variables, {result.iterations} iterations, "
    f"{result.solve_time:.0f} s, {result.status}"
)

This is Ipopt version 3.14.19, running with linear solver spral.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:        2
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        2
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.3493641e+04 0.00e+00 1.00e+02   0.0 0.00e+00    -  0.00e+00 0.00e+00   0

   1  7.5800219e+04 0.00e+00 1.23e+03   1.6 4.46e+01    -  8.45e-01 3.57e-02h  1

   2  1.1909606e+04 0.00e+00 3.97e+02   1.0 4.98e-01    -  5.66e-02 1.00e+00f  1

   3  1.3129178e+03 0.00e+00 1.14e+02   0.9 2.63e-01    -  1.00e+00 1.00e+00f  1

   4  3.4916966e+02 0.00e+00 1.70e+01   0.2 9.18e-02    -  9.98e-01 1.00e+00f  1

   5  3.4358687e+02 0.00e+00 3.81e+00  -1.2 1.92e-02    -  9.98e-01 1.00e+00f  1

   6  3.0657998e+02 0.00e+00 2.98e+00  -2.9 5.82e-03    -  1.00e+00 1.00e+00f  1

   7  3.0458641e+02 0.00e+00 7.76e+00  -4.1 2.16e-02    -  1.00e+00 1.00e+00f  1

   8  2.8312478e+02 0.00e+00 2.77e+00  -5.5 9.79e-03    -  1.00e+00 1.00e+00f  1

   9  3.0711983e+02 0.00e+00 1.67e+00  -7.2 2.79e-03    -  1.00e+00 1.00e+00h  1

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  3.3865675e+02 0.00e+00 1.50e+00  -7.4 1.57e-02    -  1.00e+00 1.00e+00h  1

  11  3.3761663e+02 0.00e+00 1.50e+00  -7.4 1.21e-02    -  1.00e+00 1.00e+00f  1

  12  3.3697149e+02 0.00e+00 1.51e+00  -7.4 1.24e-02    -  1.00e+00 1.00e+00f  1

  13  3.3644324e+02 0.00e+00 1.73e+00  -7.4 1.29e-02    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 13

                                   (scaled)                 (unscaled)
Objective...............:   2.1208406009958218e+00    3.3644323779736806e+02
Dual infeasibility......:   1.7309071566027823e+00    2.7458546758329589e+02
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   4.3669318617159297e-08    6.9275583186512660e-06
Overall NLP error.......:   1.7309071566027823e+00    2.7458546758329589e+02


Number of objective function evaluations             = 14
Number of objective gradient evaluations             = 14
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 0
Total seconds in IPOPT                               = 230.585

EXIT: Solved To Acceptable Level.

calibration: 237.3 s
shooting: 2 decision variables, 13 iterations, 231 s, Solved_To_Acceptable_Level

5. After: the fit on the training and validation days, and what the two fields became¶

In [8]:

Copied!





recovered = pd.DataFrame({"assumed (YAML)": spec.nominal, "calibrated": result.factors, "true": truth["true value"]})
recovered["error %"] = 100 * (recovered["calibrated"] - recovered["true"]) / recovered["true"]
display(recovered.round(3))
display(result.metrics().round(3))
recovered = pd.DataFrame({"assumed (YAML)": spec.nominal, "calibrated": result.factors, "true": truth["true value"]})
recovered["error %"] = 100 * (recovered["calibrated"] - recovered["true"]) / recovered["true"]
display(recovered.round(3))
display(result.metrics().round(3))

	assumed (YAML)	calibrated	true	error %
constructions[CAVITYWALL:CAL].layers[1].thickness	0.10	NaN	0.10	NaN
material[ROCKWOOL:CAL].thermal_conductivity	0.04	0.052	0.06	-13.627
spaces[SPACE:001].parameters.ach	0.30	0.628	0.60	4.654

	rmse	mae
training (Trano)	1.351	1.267
training (calibrated)	0.179	0.145
validation (Trano)	1.140	1.090
validation (calibrated)	0.148	0.127

In [9]:

Copied!





x0 = study.initial_state(result)
after = simulate_fit(calibrated, training, x0)
display(plot_fit(after, "after calibration: training days (dotted: measured)"))
after_validation = simulate_fit(calibrated, validation)
before_validation = simulate_fit(model, validation)
display(plot_fit(after_validation, "after calibration: validation days (dotted: measured)"))
x0 = study.initial_state(result)
after = simulate_fit(calibrated, training, x0)
display(plot_fit(after, "after calibration: training days (dotted: measured)"))
after_validation = simulate_fit(calibrated, validation)
before_validation = simulate_fit(model, validation)
display(plot_fit(after_validation, "after calibration: validation days (dotted: measured)"))

6. Identifiability: what the data can and cannot tell apart¶

The calibrated model reproduces the air temperature and the electricity of the true building on days it has not seen, but the two fields did not land on their true values: a leakier house and a worse insulation are both a loss proportional to the indoor-outdoor temperature difference. The data of a thermostat-controlled house fix the total loss; how it splits between the walls and the infiltration is a flat valley of the objective, and the prior (the YAML) and the first steps of IPOPT decide where along it the fit stops.

What a measurement can add: a blower-door test gives the infiltration rate. With the YAML regenerated at that rate and the conductivity as the only declared parameter, the same fit recovers it.

In [10]:

Copied!





with timed("calibration with the infiltration measured"):
    calibrated_known, result_known = study.calibrate(known_infiltration=True)
print(
    f"{result_known.method}: {result_known.variables} decision variable, {result_known.iterations} iterations, "
    f"{result_known.solve_time:.0f} s, {result_known.status}"
)
known = pd.DataFrame(
    {
        "assumed (YAML)": {study.CONDUCTIVITY: spec.nominal[study.CONDUCTIVITY]},
        "calibrated": result_known.factors,
        "true": {study.CONDUCTIVITY: study.TRUTH[study.CONDUCTIVITY]},
    }
)
known["error %"] = 100 * (known["calibrated"] - known["true"]) / known["true"]
display(known.round(4))
display(result_known.metrics().round(3))
with timed("calibration with the infiltration measured"):
    calibrated_known, result_known = study.calibrate(known_infiltration=True)
print(
    f"{result_known.method}: {result_known.variables} decision variable, {result_known.iterations} iterations, "
    f"{result_known.solve_time:.0f} s, {result_known.status}"
)
known = pd.DataFrame(
    {
        "assumed (YAML)": {study.CONDUCTIVITY: spec.nominal[study.CONDUCTIVITY]},
        "calibrated": result_known.factors,
        "true": {study.CONDUCTIVITY: study.TRUTH[study.CONDUCTIVITY]},
    }
)
known["error %"] = 100 * (known["calibrated"] - known["true"]) / known["true"]
display(known.round(4))
display(result_known.metrics().round(3))

calibration with the infiltration measured: 7.9 s
shooting: 1 decision variable, 17 iterations, 171 s, Solved_To_Acceptable_Level

	assumed (YAML)	calibrated	true	error %
material[ROCKWOOL:CAL].thermal_conductivity	0.04	0.0593	0.06	-1.2141

	rmse	mae
training (Trano)	0.344	0.301
training (calibrated)	0.168	0.138
validation (Trano)	0.297	0.237
validation (calibrated)	0.149	0.128

Key takeaways¶

The YAML is the specification of the calibration too: three fields carry {value, parameter: true, min, max}, CalibrationSpec writes the flat copy Trano reads and measures which parameters of the CasADi model each field moves; the fit adjusts the fields within their bounds and the model follows.
A closed-loop model is calibrated by simulation, not collocation: its thermostat and hysteresis make the collocated problem non-smooth; the single-shooting fit (two decision variables, CVODES inside the objective, wide central differences) brings the air-temperature error from 1.35 K to 0.18 K on the training days and from 1.14 K to 0.15 K on two unseen days, with the energy use to match.
The data decide what is identifiable: the free fit lands at a conductivity of 0.052 and 0.63 air changes per hour (true: 0.060 and 0.60) - the right total loss, split along a flat valley the objective barely feels; with the infiltration measured, the conductivity comes out at 0.0593. Declare in the YAML what a measurement cannot give, fix the rest.
The remaining 0.17 K is the switching: a thermostat that flips one integration step earlier moves the hourly sample; the fitted model has the right parameters but not the exact switching times, which is what IPOPT's jumpy log in section 4 shows.