thermal — 5R1C Thermal Model

Source:

buem/thermal/model_buem.py

Purpose

Implements the ISO 13790 5R1C (five-resistor, one-capacitor) hourly thermal model. Given a building envelope, weather, and occupancy it computes annual heating and cooling demand by solving a constrained optimisation problem.

5R1C Network

Three temperature nodes are connected by thermal conductances:

  • T_air — interior air

  • T_sur — internal surfaces

  • T_m — thermal mass

External air temperature T_e (from weather) and the HVAC heat flow Q_HC (the decision variable) close the network.

Three energy-balance equations per time-step enforce the physics:

  1. Air-node balance — convective coupling to surfaces and ventilation.

  2. Surface-node balance — radiative coupling to mass and windows.

  3. Mass-node dynamics — forward-Euler integration with annual-periodic wrap-around (\(T_m(n{-}1) \to T_m(0)\)), so no initial-condition guess is needed.

Gain Distribution (ISO 13790 §C.2)

Internal and solar gains are split across nodes:

\[\begin{split}\phi_{ia} &= 0.5\,\phi_{int} \\ \phi_{st} &= (1 - f_{Am} - f_w)\,(0.5\,\phi_{int} + \phi_{sol}) \\ \phi_m &= f_{Am}\,(0.5\,\phi_{int} + \phi_{sol})\end{split}\]

where \(f_{Am} = A_m / A_{tot}\) and \(f_w\) is the window-area fraction.

Solar Gains

Plane-of-array (POA) irradiance per surface element is computed with pvlib (isotropic-sky diffuse model).

  • Window gains: \(Q_{win} = \text{POA} \times A_{win} \times g_{gl} \times (1-F_f) \times F_w\)

  • Opaque gains: \(Q_{opaque} = \alpha \times R_{se} \times U \times A \times \text{POA}\)

The \(R_{se} \times U\) factor ensures that only the ~4–7 % of absorbed solar energy that actually penetrates inward is counted; without it cooling loads are vastly over-estimated.

Solver

LP (default) — Minimises \(\sum |Q_{\text{HC}}|\) subject to the 3 × 8760 physics constraints plus comfort dead-band \(T_{lb} \le T_{air} \le T_{ub}\). Solved via CVXPY using CLARABEL (preferred) with OSQP fallback. Positive \(Q_{\text{HC}}\) = heating, negative = cooling.

MILP (experimental) — Separates heating \(Q_h \ge 0\) and cooling \(Q_c \ge 0\) with binary indicators and big-M constraints. Solved via CVXPY (CBC / GLPK_MI) or PuLP + CBC.

Key Methods

Method

Role

__init__

Initialise model, declare containers

_initEnvelop()

Parse components tree → U-values, areas, conductances H [kW/K]

_init5R1C()

Compute network parameters (A_f, A_m, C_m, H_is, H_ms, H_win); build solar-gain profiles via pvlib

_calcRadiation()

POA irradiance per element (kW/m²)

_addConstraints()

Assemble sparse 3n × 4n equality matrix

sim_model()

Build and solve LP or MILP

_readResults()

Populate detailedResults DataFrame (T_air, T_m, T_sur, Q_HC, gains)

Design Note — LP vs. Parameterisation

Earlier BuEM versions separated heating and cooling via rule-based parameterisation. The current LP formulation unifies them into a single decision variable \(Q_{\text{HC}}\) whose sign determines the mode. This removes ad-hoc heuristics and produces sparse (zero-in-dead-band) profiles that match expected HVAC behaviour.