heatpumpmodel

Single-source implementation of the steady-state (bin-style) heat-pump performance model of

A. Rogeau, R. Vieubled, M. de la Ruche, G. Girard, "A generic methodology for mapping the performance of various heat pump configurations considering part-load behavior", Energy and Buildings 2024, https://doi.org/10.1016/j.enbuild.2024.114471

numpy is the only dependency. The package holds physics only: it knows nothing about a building stock, a weather reader or a simulation pipeline. Callers pass hourly arrays in and keep their own integration layer.

It exists so that buildingmodel and building_eload share one copy of the physics instead of two that drift apart (buildingmodel issue #47).

Install

pip install heatpumpmodel
# or, from a checkout:
pip install -e ".[dev]"

Use

import numpy as np
from heatpumpmodel import (
    Emitter, HeatPumpConfig, Mode, System, Technology,
    compute_T_base, hourly_power_split, scop, size_heat_pump,
)

# hourly series over a year: outdoor temperature (°C), relative humidity (%),
# and heat demand in any single consistent unit (it cancels in every ratio).
t_out = ...
rh = ...
demand = np.clip(20.0 - t_out, 0.0, None)

cfg = HeatPumpConfig(System.A_W, Mode.M, Emitter.MT, Technology.ON_OFF)
sizing = size_heat_pump(cfg, demand, t_out, compute_T_base(t_out))
p_h, p_e, p_h_backup = hourly_power_split(cfg, sizing, demand, t_out, rh=rh)

print(scop(p_h, p_e))            # seasonal COP

seasonal_performance(cfg, sizing, demand, t_out, rh=rh) bundles {"scop", "ecr", "peak_share"} in one call.

Configuration axes

Axis	Values
System	A_A (air/air), A_W (air/water), G_W (ground/water)
Mode	M, M_SB, BA (bivalent alternative), BP (bivalent parallel)
Emitter	FH 35 °C, LT 45 °C, MT 55 °C, HT 65 °C, FAN_COIL (sink = indoor setpoint; direct expansion only)
Technology	ON_OFF, BI_COMPRESSOR, INVERTER

Sub-models (COPCurve, TwoBranchCOPCurve, DefrostModel, PartLoadModel, WeatherCompensation) are dataclasses on HeatPumpConfig and can be replaced with manufacturer-specific fits.

The A/W COP curve (behavioural break in 0.3.0)

The default A/W full-load curve is TwoBranchCOPCurve: a defrost-degraded branch (5.60 − 0.09·ΔT + 0.0005·ΔT², T_ext ≤ −3 °C) and a frost-free branch (9.302 − 0.223·ΔT + 0.0017·ΔT², T_ext ≥ +6 °C), linearly interpolated in the outdoor-air temperature between the two. Before 0.3.0 only the degraded branch existed, applied at every outdoor temperature — the collapsed form of the source's Eq. (3.8) that the authors' published code produces. A/W SCOP therefore rises — 4.4 % to 19.5 % across the emitter × technology × mode grid on a synthetic Paris-like year, most in mild weather; A/A and G/W are unaffected.

Because heating_capacity (Eq. 6) reads the same curve, it is a capacity change too: for A/W, M mode, floor heating, inverter on that year, capacity at T_so = +10 °C is +31 %, mean capacity over heating hours +22 %, hours below minimum modulation go 843 → 1402, T_biv2 11.8 → 9.7 °C and ECR 81.1 → 89.7. Q_nom/P_e_nom are unchanged in M / M+SB mode (their design point sits on the degraded branch for any T_base ≤ −3 °C) but move slightly in BA/BP, whose backup bisection runs on the branch-aware capacity.

The best case of that grid (FH, inverter, BA) now reaches SCOP ≈ 4.67, at or above the top class of the paper's own European map — headroom that follows from implementing the printed source faithfully, not from any empirical validation. To reproduce pre-0.3.0 A/W numbers (and to reproduce the paper's Fig. 4, which is an output of the collapsed curve):

from heatpumpmodel import HeatPumpConfig, LEGACY_A_W_SINGLE_CURVE, System, Mode, Emitter, Technology

cfg = HeatPumpConfig(System.A_W, Mode.M, Emitter.MT, Technology.ON_OFF,
                     cop_curve=LEGACY_A_W_SINGLE_CURVE)

TwoBranchCOPCurve.cop_fl requires the source (outdoor-air) temperature and raises if it is missing — so COP_CURVES[System.A_W].cop_fl(dT) without t_so, which used to return an array, now raises. See its docstring in heatpumpmodel/core.py for the provenance chain (PhD thesis tel-02969503, not Ruhnau et al.) and for the documented overlap between the branch structure and DefrostModel — which also means an A/W-vs-A/A SCOP gap now straddles two different defrost representations and partly reflects curve lineage, not physics.

Check the sizing solve status on BA/BP

sizing = size_heat_pump(cfg, demand, t_out, compute_T_base(t_out))
if not sizing.converged:                 # sizing.solve_status names which rail was hit
    ...                                  # surface it; do not publish a railed Q_nom

For Mode.BA / Mode.BP, Q_nom is solved by bisection for a 10 % backup share. When that target is unreachable the solver returns a bracket rail (4·h_base or 1e-3·h_base) as a best-effort size, and every quantity derived from it is an artefact of the search interval rather than a sizing. Since 0.3.0 HeatPumpSizing.solve_status (SizingSolveStatus.CONVERGED / RAILED_LOW / RAILED_HIGH) and the derived converged property make that visible.

Since 0.4.0 a fourth status, MAX_ITER, covers the bisection's third silent exit: the target was bracketed but |achieved − target| never met the tolerance, and the final midpoint is returned — the exit that looks most like a normal answer (interior value, no rail, no exception). converged is False there. It is a reporting change only: on the committed 13 824-solve grid (scripts/differential_sizing_status.py, which regenerates these numbers), 756 solves (all BA) reclassify and no number moves. The flip count and the residual ceiling are properties of that grid — a different but equally defensible grid moves both; what is grid-independent is the 0 numeric moves, the confinement to BA, and residuals of the order of the solver's tolerance. Consumers wanting the old grouping test solve_status in (CONVERGED, MAX_ITER). converged still means only "no rail and no max_iter exhaustion", and it is also the default for the closed-form M / M+SB rule, so check Q_nom for plausibility as well.

Declared capacity, and the backup-share target

from heatpumpmodel import DeclaredCapacity, SizingBasis, compute_T_base, size_heat_pump

sizing = size_heat_pump(cfg, demand, t_out, compute_T_base(t_out),  # 8 kW at 55/7 °C
                        declared=DeclaredCapacity(8000.0, sink_c=55.0, source_c=7.0))
sizing.basis is SizingBasis.DECLARED     # T_biv / T_biv2 are NaN — see below

declared skips the sizing rule and uses the caller's capacity, rescaling P_e_nom and P_e_min through the package's own Eq. 5 / Eq. 13 at the declared rating point. The rating point is required, not defaulted: a declared thermal power is not automatically Q_nom (G/W nameplates are rated at 0 °C source water — a different physical quantity from an air-source rating). The bivalence diagnostics are properties of the derived solve, so they are invalidated (NaN) rather than carried, and the choice is recorded on sizing.basis (DERIVED by default; every existing call is unchanged). Neither convention is endorsed here — on the ADEME campaign derived and declared arms err in opposite directions on backup engagement and neither reproduces the metered value.

size_heat_pump(..., backup_target=...) promotes the solver's private default to the public signature (validated to (0, 1), fail-loud outside it). It defaults to None, the sentinel for "no target was asked for", which takes the rule's own 0.10 convention — unchanged and bit-identical. The 10 % share is a convention of the sizing rule, not a prediction, and it has to be reachable by a caller comparing against measured data. Passing declared together with any explicit backup_target raises, 0.10 included: the caller asked for a target and the declared path could only ignore it.

Meter-boundary layers: auxiliaries and DHW (0.4.0)

heatpumpmodel.core predicts machine electricity. What a dwelling's distribution board draws is more than that (heatpumpmodel#5 §1.3):

E_total = E_compressor(COP*) + E_backup            <- core
        + E_aux(standby, circulators, source pump) <- heatpumpmodel.auxiliaries
        + E_DHW(charge conversion + boost)         <- heatpumpmodel.dhw

Both are additive layers: they report their own terms beside the core's, never folded into them, and they move nothing the core computes — no space-heating COP, no capacity, no sizing. A consumer that never calls them gets the same numbers byte for byte.

auxiliaries predicts the electricity inside a declared metering perimeter, and the perimeter is a derived output, not a label the caller asserts: plan_aux_terms turns three recorded installation facts (ground source? separate auxiliaries circuit? decoupling bottle or buffer tank?) into the inclusion set and an H4* / H4*-ext / H1*+bu label. A consumer wired outside the perimeter is suppressed, not added — an "auxiliaries" breaker reduces the prediction, because the meter being predicted cannot see what is on it. Results come back component-wise (standby, circulator, overlap deduction, source pump) for three named literature coefficient cases; there is no default coefficient set.

dhw converts a DHW charge-heat series into the electricity the same compressor spends making it: the same machine at a different operating point (scope A). DHW demand, draw profiles and tank dynamics are scope B and are not modelled. The sink is the charge temperature, passed explicitly, so the space-heating sink and the sizing solve cannot move; the immersion boost is a measured pass-through, never a prediction.

Status, in the words that must travel with any number from these layers: the auxiliary circulator coefficient is an open question — it runs +15 … +51 % above its own campaign's same-sample energy anchor (criterion K13 FAIL), and it is reported, never tuned. That comparison is measured on the campaign's mild year: R8 states a 2021-like winter would cut SCOP by 15–20 %, and an auxiliary share of electricity is mechanically higher in a mild year on both sides of an anchor comparison — so the vintage does not explain the overshoot away, and equally nothing here generalises to a design winter. The DHW layer is prototyped and gated as a diagnostic, not validated (the pre-registered rule returns FAIL on the total clause), quotable only with its target's own 1.00–1.13× one-directional fidelity, measured on 4 dwellings, in the same breath — and those four are Δ ranks 24, 27, 31 and 34 of 34, i.e. an agreement on the high-Δ tail rather than on a random four. On the pre-registered basis the informative subset is 1 pass (the strawman clause) / 5 fail, and on the split basis 2 pass / 1 fail / 3 uninformative; it is never "passes five of six". Read the two module docstrings (heatpumpmodel/auxiliaries.py, heatpumpmodel/dhw.py) before publishing either: they carry the provenance, the guards, the counted-flag semantics and the full status wording.

From French DPE data

heatpumpmodel.dpe maps buildingdata's heat_pump_type / heating_emitter_type / heat_pump_installation_period DPE columns onto a HeatPumpConfig, so buildingmodel and building_eload resolve the same building to the same machine in both the static and dynamic stages instead of each guessing independently:

from heatpumpmodel import HeatPumpConfig, System, Mode, Emitter, Technology, config_from_dpe

default = HeatPumpConfig(System.A_W, Mode.M, Emitter.MT, Technology.ON_OFF)
cfg = config_from_dpe("air/air", "air", "[2015, 2100]", default=default)

Any DPE column that is None falls back to the matching axis of default; a present value outside the module's vocabulary raises ValueError rather than drifting silently. See heatpumpmodel/dpe.py for the vocabulary vintage and the full mapping tables.

The emitter is resolved against the system, not from heating_emitter_type alone. The DPE's "air" emission class is two different machines: on an air/air generator it is direct expansion and the sink is the room (FAN_COIL), while on a hydronic generator it is a fan-coil water loop — a ventilo-convecteur, still circulating water at roughly 35–45 °C — which resolves to LT (45 °C). The invariant is FAN_COIL ⇔ A_A. Before 0.5.0 "air" mapped to FAN_COIL unconditionally, sinking hydronic units to the indoor setpoint and overstating their SCOP (up to 10.3 on ground-source records); that fix moves the numbers for ~15 % of the DPE heat-pump stock, so a consumer upgrading from 0.4.x must re-measure rather than bump the floor silently (heatpumpmodel#6).

The same invariant is enforced at the constructor: since 0.5.0 HeatPumpConfig(System.A_W | System.G_W, ..., Emitter.FAN_COIL, ...) raises ValueError instead of silently sinking the machine to the indoor setpoint. Reach for Emitter.LT for a hydronic fan-coil loop, or System.A_A if the unit really is direct expansion.

Conventions

• Temperatures in °C, ΔT gaps in K.
• Powers/demand in one arbitrary, self-cancelling unit.
• Relative humidity in percent [0, 100] — the EPW convention.
• Air-source configs (A/A, A/W) require an rh series: a missing, all-NaN, all-zero or partially non-finite rh, a [0, 1] fraction-convention series, or any value outside [0, 100] raises ValueError rather than silently skipping the defrost derate (which would leave SCOP ~5 % optimistic). The checks are whole-series: one bad summer hour fails the call.
• Every effective COP is floored at 1.0 — a heat pump never draws more electricity than the resistance backup would for the same heat. The floor is exported as HP_COP_FLOOR, so layers and consumers test against it instead of mirroring the literal.

Documentation

doc/heat_pump_model_spec.md is the implementation contract: equation-by-equation mapping to the paper, coefficient provenance (including the values resolved from the authors' Zenodo code rather than the PDF), and the documented deviations.

Tests

pytest                          # hermetic suite, synthetic climate
HEATPUMPMODEL_PARIS_EPW=/path/to/paris.epw pytest -m integration

The integration test reproduces the paper's Fig. 4 SCOP values and needs the authors' Paris-Montsouris TMY EPW; it skips when that file is not supplied. See its docstring for why an ERA5 Paris record does not substitute.

Licence

MIT — see LICENSE.