Adjusted layered PTES model

[cpschau]

Layered-PTES model notes — formulation, wiring, and the comparison with Amos' model

This documents the current, energy-exact layered-PTES model (the toy's
recharge model, now also ported into PyPSA-Eur:
scripts/prepare_sector_network.py, scripts/build_ptes_operations/,
scripts/solve_network.py) and compares it, component by component, with the
Amos model from PyPSA-Eur PR
#2129 ("feat: layered ptes
model"), which the port replaces.

The two models share the same idea — a stack of fixed-temperature water layers
(one PyPSA Store per layer, volume in m³, energy density
m3_to_mwh[l] = ρc_p(T_l − T_bottom)) charged from and discharged to the
urban-central heat bus, with a booster heat pump for sub-forward layers. They
differ in how volume and energy are tracked, and that difference is the whole
story: the current model conserves both exactly, Amos's creates energy on
discharge.

---

1. Preliminaries — energy reference, temperatures, ΔT and COP

1a. Notation and the "bottom-referenced store debit"

symbol	meaning
T_l	fixed temperature of layer l (exogenous, hottest-first)
T_bottom	coldest (bottom) layer temperature
T_ret, T_fwd	DH return / forward temperature (per node, time-varying)
T_retL	temperature of the discrete return-level layer (the layer chosen to represent the return)
T_boost	temperature of the discrete reinjection / boost-deposit layer (where boosted, cooled water is put back; = layer hp_return_layer)
e_l	state of charge of layer l's Store in m³
Δ(a,b)	ρc_p(T_a − T_b) — enthalpy per m³ between two temperatures [MWh/m³]

A layer Store holds volume (m³), not energy. Each layer is given a fixed
energy density measured relative to the bottom layer,

m3_to_mwh[l] = ρc_p (T_l − T_bottom)         (so m3_to_mwh[bottom] = 0)

so the energy content of the whole PTES is the bottom-referenced sum
E = Σ_l m3_to_mwh[l]·e_l. The choice of reference is arbitrary because every
operation only moves volume between layers: moving 1 m³ from layer a to
layer b changes E by m3_to_mwh[a] − m3_to_mwh[b] = ρc_p(T_a − T_b) = Δ(a,b),
in which the bottom term cancels. We call that change the store debit (when
E falls, on discharge) or store credit (when E rises, on charge). The
whole point of the current model is that, on every operation,

heat drawn from DH (charge)        =  store credit                          (charge)
heat delivered to DH − electricity =  store debit  =  Σ Δ(from, to) per m³   (discharge)

i.e. delivered useful energy is backed exactly by a bottom-referenced volume
relocation, with no leftover.

1b. The heat-pump output relation (shared, and correct in both models)

A discharged m³ of layer water at T_src gives up heat as it cools. Part is
delivered directly (heating the DH return up to T_src); the sub-return part
is the evaporator heat the booster lifts to T_fwd:

q_hp_out   = COP/(COP−1) · q_evaporator                   (booster output)
q_released = q_direct + q_evaporator = Δ(src, deposit)·v  (enthalpy released)
⇒  q_delivered_to_DH − q_elec = q_released                (energy balance)

Both models get this relation right; they differ in how the cooling depth ΔT and the COP are
parametrised (§1c) and in how the volume is booked (§2–§4).

1c. Per-layer cooling depth ΔT, reinjection layer, and COP

A booster discharging layer n cools the water from its source-inlet temperature
down to a deposit temperature, and that cooling depth fixes both where the
spent volume is reinjected and how efficient the heat pump is.

• Cooling depth / reinjection target. The ideal evaporator outlet is
  T_target = T_ret − ΔT_HP with ΔT_HP = (T_fwd − T_source)·(1 − 1/COP)
  (T_source = min(T_n, T_ret)), i.e. the booster cools the return-level water
  below the return temperature by the amount the COP can repay. This continuous
  target is snapped to a discrete layer to give hp_return_layer ( = the
  T_boost layer):

  • Amos: the nearest layer to T_target (strictly colder than the source).
  • Current model: with conservative_return_layer, the warmest layer at or
    below T_target (a floor), so the deposit is never warmer than intended
    and discharge cannot create energy; otherwise nearest, with a warning when a
    layer ends up above the return temperature.

• COP per layer per timestep. The Jensen et-al. thermodynamic model
  (CentralHeatingCopApproximator) is evaluated with sink_outlet = T_fwd,
  sink_inlet = max(T_n, T_ret) (after preheating), source_inlet = min(T_n, T_ret)
  and a source_outlet:

  • Amos: source_inlet − heat_source_cooling (a fixed shallow ΔT, e.g. 6 K)
    ⟹ optimistic, the COP ignores how deep the water is actually cooled.
  • Current model: source_outlet = T[hp_return_layer] (the real evaporator
    depth) — recomputed in PtesApproximator.booster_cop and exported as
    booster_cop. Deeper cooling ⟹ lower COP ⟹ more booster electricity, honestly.
    Energy conservation does not depend on the COP value (heat = evaporator +
    elec for any COP); the deposit-depth COP only sharpens the elec/source split.

• The ΔT–COP circularity is broken, not solved. Note ΔT_HP depends on the COP
  and the COP (at source_outlet = T_ret − ΔT_HP) depends on ΔT_HP — an implicit
  equation. It is not solved as a fixed point. Instead a one-shot, two-pass
  forward scheme is used: (1) build_cop_profiles computes COP₀ once with a
  fixed shallow source cooling, so COP₀ carries no ΔT_HP dependence and the loop
  is cut; (2) hp_return_layer evaluates ΔT_HP = (T_fwd−T_source)(1−1/COP₀) per
  timestep, takes the min over time, and snaps it to a discrete layer; (3)
  booster_cop re-evaluates COP₁ once at that snapped depth, and that is the COP
  the LP uses. The COP that picks the layer (COP₀) and the COP the LP runs
  with (COP₁) are deliberately left slightly inconsistent — it does not matter,
  because the deposit is discrete anyway and (per §1a) the LP energy balance is
  COP-independent, so the residual only nudges the electricity/source split, never
  the energy total. A self-consistent value would be a scalar Picard iteration
  ΔT ← (T_fwd−T_source)(1−1/COP(T_ret−ΔT)) per (layer, node), but it is not worth
  the cost.

---

2. Wiring of the current model (energy-exact, volume-conserving)

Per layer n, the PyPSA components and their link efficiencies (notation, incl.
Δ(a,b), T_retL, T_boost, in §1a). nb = needs_boosting (1 when
T_fwd > T_n, else 0). The three cooling depths are clipped at zero:

direct_drop = max(T_n − T_retL, 0)        above-return part, delivered directly
below_drop  = max(T_retL − T_boost, 0)    sub-return part, lifted by the booster
extract_dt  = direct_drop + below_drop = max(T_n − T_boost, 0)   total cooling

2a. Charging — volume trade (no volume created)

                         ┌──────────────────────────────┐
   URBAN CENTRAL HEAT ───┤ water-pit charger (one per    │
        bus       q ────►│ cold source s, hot dest d)    │
     (charger bus0)      │                               │
                         │  bus1 → Store water pits d    │  eff  = 1 / Δ(T_d,T_s)
                         │  bus2 → Store water pits s     │  eff2 = −eff
                         └───────────────────────────────┘
   one link per (source layer s with T_s ≤ T_ret, dest layer d with T_d > T_ret).
   p_max_pu = charging_availability(d) · 1{ s ≥ return-layer > d }   (per node)

+x m³ appear at d, −x m³ leave s ⇒ volume conserved. Heat drawn
= Δ(T_d,T_s)·x = stored-energy gain m3_to_mwh[d]−m3_to_mwh[s] ⇒ exact.

2b. Discharging — discharger + the direct/boost split

   Store water pits n
        │  discharger      eff  = 1  → bus1  (+1 m³ to the RETURN-level layer Store)
        │  (volume)        eff2 = 1  → bus2  (+1 token to the preheater bus)
        ├───────────────────────────────► Store water pits  RETURN layer
        └───────────────────────────────► ptes preheater n  (m³ token bus)
                                                 │
              ┌──────────────────────────────────┴───────────────────────────────┐
              │ DIRECT-UTIL  (no boost)                 │ BOOST multilink          │
              │ p_max_pu = 1 − nb                        │ p_max_pu = nb            │
              │ bus0 = token                             │ bus0 = token             │
              │ bus1 = heat   eff = Δ(T_n,T_returnLayer) │ bus1 = RETURN layer eff =−1 (−1 m³)
              ▼                                          │ bus2 = BOOST  layer eff2=+1 (+1 m³)
      URBAN CENTRAL HEAT                                 │ bus3 = hp source eff3 = Δ(T_n,T_boost)
                                                         ▼
                                              ptes hp source n  (MWh bus)
                                                         │  HEAT-EXCHANGER split
                                              ┌──────────┴──────────┐
                                  direct_frac │                     │ hp_frac
                              = direct_drop   │                     │ = below_drop
                                / extract_dt  │                     │   / extract_dt
                                              ▼                     ▼
                                     URBAN CENTRAL HEAT     ptes hp inlet n  (MWh bus)
                                                                    │  BOOSTER heat pump (p ≤ 0)
                                                                    │  bus0 = heat        (q_out = −p0)
                                                                    │  bus1 = electricity eff  = 1/COP
                                                                    │  bus2 = hp inlet    eff2 = (COP−1)/COP
                                                                    ▼
                                                       URBAN CENTRAL HEAT   (+ draws elec from
                                                                             the low-voltage bus)

The split fractions always sum to 1 (when extract_dt > 0; set to 0 otherwise to
avoid 0/0), and their boundary cases say which path the discharge takes:

• direct_frac = 1, hp_frac = 0 ⟺ below_drop = 0 ⟺ T_retL = T_boost: there
  is no usable boosting depth — the return-level layer is already the deposit
  layer (common under the conservative floor, where T_retL is pushed to the
  bottom). The booster never runs and the whole release goes direct. This is why
  the solved 5-layer case draws ~0 booster electricity.

• direct_frac = 0, hp_frac = 1 ⟺ direct_drop = 0 ⟺ T_n ≤ T_retL: a
  sub-return layer is being cooled, so the entire extraction is evaporator heat
  (no above-return part) — see case (iv).

• 0 < direct_frac < 1: a genuinely boosted warm layer (T_n > T_retL > T_boost),
  with both a direct part and an evaporator part.

• NB on bus order: electricity must sit on bus1 of the booster — PyPSA-Eur's
  distribution-grid rewiring keys off bus1 for heat-pump electricity and would
  otherwise append " low voltage" to whatever is on bus1, mangling the
  evaporator bus into a phantom free-energy source.

2c. Stores, interlayer links, generators, constraints

   Store water pits l  (m³, e_cyclic, standing_loss = 0)            for each layer l
        │  interlayer link  bus0=l → bus1=l+1   (downward volume flux = standing loss)
        ▼
   Store water pits l+1
   ──────────────────────────────────────────────────────────────────────────────
   Store water pits  (aggregate, e_max_pu, carries the €/m³ storage capital_cost)
   Link  ptes heat pump  (dummy, p_max_pu=0, carries the booster €/MW capital_cost)
   GENERATORS:  none on the PTES buses  (no vents, no bottom-water source)

   EXTRA CONSTRAINTS (solve_network.py):
     volume capacity      Σ_l e_l(t) = ē / m3_to_mwh[top]            (∀t; tank always full)
     interlayer flow      p_{l→l+1}(t) = Ψ_l · e_l(t)                (standing-loss channel)
     HP capacity          Σ_l p_nom(booster_l) ≤ p_nom(dummy HP)     (shared investment)
     throughput           Σ chargers + Σ m3_to_mwh·dischargers ≤ R·ē

Everything references the discrete deposit layers, so on every discharge
delivered − elec = ρc_p(T_n − T_deposit) = the bottom-referenced store debit.
Volume is conserved at every link. The only dissipation is the interlayer flux
(standing_loss = 0 on the stores, so the loss is single-channel).

---

3. Wiring of the Amos model (PR #2129) — where energy leaks in

Same layer stores, but charging mints volume and discharging is a
preheater → heat-pump chain whose heat references the actual return temperature
while its volume is booked against discrete layers. Same Δ/nb notation as §2.

3a. Charging — create-volume (volume minted, not traded)

                         ┌──────────────────────────────┐
   URBAN CENTRAL HEAT ───┤ create-volume charger         │  eff = charging_avail[n] / m3_to_mwh[n]
        bus       q ────►│ (one per layer n)             ├──► Store water pits n   (+ q·eff m³,
     (charger bus0)      │                               │                          CREATED from heat)
                         └───────────────────────────────┘
   1 MWh of heat ⟹ 1/m3_to_mwh[n] m³ minted at layer n (stored energy = 1 MWh,
   bottom-referenced). No source layer is drained ⇒ Σ_l e_l is NOT conserved here;
   it is rebalanced only by the free generators below.

3b. Discharging — discharger → preheater → heat pump

   Store water pits n
        │  discharger   eff = 1   (moves +1 m³ as a unit onto the preheater bus)
        ▼
   ptes preheater n   (m³ bus)
        │   PREHEATER multilink   (consumes the 1 m³ unit)
        ├── bus1 = heat                    eff  = Δ(T_n, T_ret)        ◄ delivers to DH at the
        │                                                                ACTUAL return T
        ├── bus2 = hp input                eff2 = nb                   ──► +nb   m³ to the booster
        └── bus3 = preheater_return_layer  eff3 = 1 − nb               ──► +(1−nb) m³ to a
                          │                  │                              DISCRETE layer Store
                          ▼                  ▼
                 URBAN CENTRAL HEAT   Store water pits  preheater_return_layer

   ptes hp input n   (m³ bus)
        │   HEAT PUMP   (p ≤ 0)
        ├── bus0 = heat                     q_out = heat_pump_eff = Δ(T_fwd, T_n)   ◄ the LIFT,
        │                                                                            not the cooling
        ├── bus1 = electricity              eff  = 1/COP
        ├── bus2 = hp input                 eff2 = 1/heat_pump_eff
        └── bus3 = hp_return_layer          eff3 = −1/heat_pump_eff    ──► +1 m³ to a DISCRETE layer
                          │                  │                  │
                          ▼                  ▼ (draws elec)     ▼
                 URBAN CENTRAL HEAT     low-voltage bus   Store water pits  hp_return_layer

3c. Stores, generators, constraints

   Store water pits l  (m³, e_cyclic, standing_loss > 0)            for each layer l
        │  interlayer link  bus0=l → bus1=l+1   (constraint p = Ψ·e, FLAT Ψ = 0.001)
        ▼
   Store water pits l+1
   ──────────────────────────────────────────────────────────────────────────────
   Store water pits  (aggregate, e_max_pu, €/m³ capital_cost);  dummy ptes heat pump (€/MW)
   GENERATORS on the layer stores:
        ptes vent (per layer)     Store l        ◄── p ≤ 0, marginal_cost = −1  (destroys volume,
                                                                                  rewarded, depth-graded)
        bottom-water generator    Store bottom   ◄── p ≥ 0, marginal_cost = +1  (creates volume)
   EXTRA CONSTRAINTS: volume capacity; interlayer flow (flat Ψ); HP capacity (Σ p_nom — sum-of-peaks).
   Stores ALSO carry a PyPSA standing_loss ⇒ standing loss is double-counted (interlayer + store).

Two structural mismatches make this leak (cf. the current model, which has neither):

1. Volume is minted/destroyed for free. The charger mints volume from heat; the
   per-layer vents (rewarded, marginal_cost = −1) destroy it; the bottom-water
   generator mints more. Σ_l e_l is not conserved per operation, so the optimiser
   can route around the physics.
2. Heat references the actual T_ret, volume references a discrete layer. The
   preheater delivers Δ(T_n, T_ret_actual), but the spent volume is debited only
   down to the nearest preheater_return_layer / hp_return_layer. When that
   discrete layer sits above the actual return temperature the discharge debits
   too little ⇒ every m³ conjures ρc_p(T_layer − T_ret_actual); and the HP delivers
   the lift Δ(T_fwd, T_n) regardless of how deep the water was actually cooled,
   so the lift and the volume debit only reconcile if hp_return_layer snaps exactly
   to T_ret − ΔT_HP.

---

4. The four cases — volume & energy, current vs Amos

Notation in §1a (T_retL, T_boost, Δ); v = 1 m³ discharged. "Store debit"
is the bottom-referenced energy drop Σ Δ(from,to)·v (§1a).

(i) Charging

	volume	energy
current	trade: −v at a sub-return source, +v at the hot dest → conserved	heat in = ρc_p(T_d−T_s)·v = stored gain → exact
Amos	charger creates +v at the dest from nothing; balanced later by free vents / bottom-water	1 MWh heat → 1 MWh stored (bottom-ref); but the non-physical volume balance is what later leaks

(ii) Direct discharge, no boost (warm layer, T_fwd ≤ T_n)

	volume	energy
current	−v at n, +v at the return-level layer → conserved	delivered ρc_p(T_n−T_retL) = debit m3_to_mwh[n]−m3_to_mwh[retL] → exact
Amos	−v at n, +v at preheater_return_layer (discrete)	delivered ρc_p(T_n−T_ret_actual), debit ρc_p(T_n−T_prl); mismatch ρc_p(T_prl−T_ret) → free lunch if T_prl > T_ret

(iii) Boosted discharge from a layer with T_n > T_ret (T_fwd > T_n)

The discharged m³ is cooled in two stages: directly down to the return level
(Δ(T_n,T_retL) straight to DH), then by the booster from there down to the
reinjection layer (evaporator Δ(T_retL,T_boost), lifted to DH at COP/(COP−1)).

	volume	energy
current	−v at n, +v at the boost layer (return level is a transient pass-through) → conserved	delivered − elec = Δ(T_n,T_boost) = debit m3_to_mwh[n]−m3_to_mwh[boost] → exact (the hp-source bus carries the full released enthalpy, the HX splits it)
Amos	−v at n, +v at hp_return_layer (discrete)	delivered = Δ(T_n,T_ret)+Δ(T_fwd,T_n); exact only if T_hp_return = T_ret − (T_fwd−T_n)(1−1/COP) lands exactly on a layer; the discrete snap leaks

(iv) Discharge of a sub-return layer (T_n ≤ T_ret)

	behaviour
current	such layers are the cold reservoir: charging draws volume up from them, discharge deposits cooled volume into them. They deliver no direct heat (Δ(T_n,T_retL) = 0 once floored); if cooled further they are pure HP-evaporator source. Volume and (low, near-bottom) energy stay exactly tracked.
Amos	preheater_efficiency = 0 for T_n ≤ T_ret, so no direct heat; the volume still flows through the create/vent/bottom-water machinery, so the cold end is where most of the free volume is minted and vented.

---

5. Net energy balance (legacy intermediate model omitted)

Toy, single node, identical boundary conditions (created = load + losses − elec; SCOP = load/elec):

model      objective   elec    load   interlayer   created   SCOP
current    32.99 M€    253.6   245.1     6.3         −2.3     0.97   (≈ conservative)
Amos       24.33 M€    167.3   245.1     7.9        +85.7     1.47   (energy conjured)

Workflow, PyPSA-Eur 5-layer scenario, 26 urban-central nodes (created = net PTES heat to DH − elec drawn + interlayer loss):

model                       created [GWh/yr]   volume drift   objective
current (this port)              0.00            5e-7         5.145e11   ← exact
Amos-style, nearest layers   +294,000            (n/a)        5.08e11    ← large free lunch

The current model nets to exactly zero energy created with volume conserved
to 5e-7; its only loss is the interlayer standing-loss channel. Amos's structure
conjures energy on discharge, which the optimiser then exploits (the objective is
artificially low and the PTES is over-cycled).

---

6. Discretisation bias (both models)

The model is a piecewise-constant approximation of a continuously stratified PTES.
For evenly-spaced layers the nominal layer temperature is the band mean, and
the Jensen COP uses both inlet and outlet temperatures (the glide) — already a
mean-temperature representation, more accurate than a single-mean-T COP. The
residual error is quantisation (collapsing each band to a point), bounded by ≈ ½ a
layer spacing. In the toy the n_layers sweep gives 3→33.5, 5→32.7, 9→31.4 M€:
coarser layers under-value the PTES (higher cost), converging downward as N
grows — i.e. conservative. Room to improve within linearity: more / non-uniform
layers concentrated near the forward and return temperatures (layer temps must
stay exogenous — endogenous temperatures are nonlinear).

Conservative discretisation (current model). Because the current model
references discrete deposit layers, it never creates energy regardless of layer
placement; the only question is whether the direct part over- or under-delivers
relative to the physical return temperature. The optional
sector.district_heating.ptes.conservative_return_layer flag floors the
return-level and reinjection layers to the warmest layer at/below the target, so
discharge stays on the conservative side; build_ptes_operations additionally
warns for any node whose return-level layer sits above its return temperature.

---

7. Implementation notes (current model)

• Charge = ~L²/2 volume-trade charger links per node (one per
  cold-source × hot-dest pair); the single-tank num_layers == 1 case keeps the
  legacy create-volume charger.
• Booster COP is recomputed at source_outlet = T[hp_return_layer] (the real
  evaporator depth) in PtesApproximator.booster_cop and exported as
  booster_cop in ptes_operations; prepare_sector_network uses it for the
  booster. Energy conservation is independent of the COP value (heat = evaporator
  • elec for any COP); this only sharpens the electricity/source split.
• Standing loss is single-channel: layer-store standing_loss = 0, the loss
  lives entirely in the interlayer-flow constraint. A nonzero store standing_loss
  would, with volume-conserving charging and e_cyclic, drive Σ_l e_l → 0.
• Booster bus order: electricity on bus1, evaporator (hp-inlet) on bus2
  (see §2b note).