Rung 39: a negative relative growth adds nothing to a carrier's growth limit — PyPSA clips it at zero¶

description: A plain `n.optimize()`, and its multi-period and stochastic classes, in one file. Every second-stage
  quantity spans a `scenario` (a future dispatch is chosen in) and every asset stands in the investment
  `period`s its build year and lifetime span. A parameter spans `scenario` exactly when PyPSA reads it
  per scenario. Capacity is chosen once, before the future is known, and paid once per active period at
  its cost in expectation over the scenarios; operation is the expectation over the scenarios' weights,
  with a share priced at the tail through the CVaR rows, which stand only where that share is positive.
  A plain run feeds one scenario, one period, all-active masks and unit weights, and the model collapses
  to the standard one. A security-constrained run copies each branch flow limit once per outage in an
  `outage` set that a plain run leaves empty. Which snapshots an asset is active in, a scenario's weight,
  and the outage factors are data prep.
dimensions:
  scenario: {description: 'the futures dispatch is chosen in, each with a weight'}
  snapshot: {description: dispatch periods, dtype: datetime}
  bus: {description: network nodes}
  generator: {description: 'generating units, each on one bus'}
  load: {description: 'demands, each on one bus'}
  store: {description: 'pure energy stores, each on one bus'}
  period: {description: investment periods — PyPSA's `investment_periods`, dtype: int}
  carrier: {description: 'energy carriers, what a growth limit is set per'}
relations:
  snapshot_period: {description: the investment period a snapshot falls in, key: snapshot, values: period}
  Generator_bus: {description: the bus a generator sits on, key: generator, values: bus}
  Load_bus: {description: the bus a load sits on, key: load, values: bus}
  Store_carrier: {description: the carrier a store holds, key: store, values: carrier}
  Store_bus: {description: the bus a store sits on, key: store, values: bus}
parameters:
  snapshot_weightings_objective:
    description: PyPSA's `snapshot_weightings.objective` — hours a snapshot stands for in the cost
    dims: [snapshot]
  Generator_p_nom:
    description: nominal power
    dims: [scenario, generator]
  Generator_p_nom_extendable:
    description: whether the nominal power is a decision
    dims: [generator]
    dtype: bool
  Generator_p_min_pu:
    description: least output, per unit of nominal power
    dims: [scenario, snapshot, generator]
  Generator_p_max_pu:
    description: most output, per unit of nominal power — an availability profile
    dims: [scenario, snapshot, generator]
  Generator_marginal_cost:
    description: cost of one unit of output
    dims: [scenario, snapshot, generator]
  Generator_marginal_cost_quadratic:
    description: cost of the square of one unit of output
    dims: [scenario, snapshot, generator]
  Generator_sign:
    description: the sign output enters its bus's balance with — PyPSA's `sign`, `1` unless given, `-1`
      for a unit that draws power. PyPSA refuses one that differs by scenario (`consistency.py:1187`)
    dims: [generator]
  Generator_committable:
    description: whether output is gated by an on/off status decision
    dims: [generator]
    dtype: bool
  Load_p_set:
    description: demand
    dims: [scenario, snapshot, load]
  Load_sign:
    description: the sign a load's demand enters its bus's balance with — PyPSA's `sign`, `-1` unless
      given, `1` for a load that feeds its bus. PyPSA refuses one that differs by scenario (`consistency.py:1187`)
    dims: [load]
  Load_active:
    description: whether a load stands in the model — PyPSA's `active`. A load has no build year and no
      lifetime, so the flag holds in every snapshot. PyPSA refuses one that differs by scenario (`consistency.py:1195`)
    dims: [load]
    dtype: bool
  scenario_weight:
    description: PyPSA's `scenario_weightings.weight` — the probability of a future
    dims: [scenario]
  CVaR_omega:
    description: PyPSA's `risk_preference['omega']` — the share of operating cost priced at the tail rather
      than in expectation; zero recovers the risk-neutral model
    dims: []
  period_weight_objective:
    description: PyPSA's `investment_period_weightings.objective` — what a period's cost weighs
    dims: [period]
  Generator_active:
    description: whether a generator stands in a snapshot's period — PyPSA's `active`, from build year
      and lifetime, data prep
    dims: [snapshot, generator]
    dtype: bool
  Store_active:
    description: whether a store stands in a snapshot's period — PyPSA's `active`, data prep
    dims: [snapshot, store]
    dtype: bool
  Store_capital_weight:
    description: the sum of period weights a store stands in — PyPSA's `active * period_weighting`, summed,
      data prep
    dims: [store]
  Store_first_active:
    description: one in the first period a store stands in, zero elsewhere, data prep. PyPSA `1.3.0` takes
      `active.cumsum() == 1`, which also counts a store that has retired in every later period (`global_constraints.py:276`,
      PyPSA/PyPSA#1938)
    dims: [period, store]
  Carrier_max_growth:
    description: most capacity of a carrier that may be added in a period; no value means no limit. The
      least over the scenarios, as PyPSA takes it (`global_constraints.py:226-230`), data prep. PyPSA
      reads it only under `multi_investment_periods` (`global_constraints.py:219-220`), so data prep feeds
      no value otherwise
    dims: [carrier]
  Carrier_max_relative_growth:
    description: share of the previous period's additions that may be added on top — the least over the
      scenarios, as PyPSA takes it, data prep
    dims: [carrier]
  snapshot_weightings_stores:
    description: PyPSA's `snapshot_weightings.stores` — hours a snapshot stands for in a storage balance
    dims: [snapshot]
  Store_e_nom_min:
    description: least nominal capacity an extendable store may be built at
    dims: [scenario, store]
  Store_e_nom_max:
    description: most nominal capacity an extendable store may be built at
    dims: [scenario, store]
  Store_capital_cost:
    description: cost of one unit of nominal capacity — PyPSA's `capital_cost`, periodized as an annuity
      in data prep
    dims: [scenario, store]
  Store_e_nom_extendable:
    description: whether the nominal energy capacity is a decision
    dims: [store]
    dtype: bool
  Store_e_min_pu:
    description: least energy held, per unit of nominal capacity — negative for a store that may go short
    dims: [scenario, snapshot, store]
  Store_e_max_pu:
    description: most energy held, per unit of nominal capacity
    dims: [scenario, snapshot, store]
  Store_sign:
    description: the sign the power a store delivers enters its bus's balance with — PyPSA's `sign`, `1`
      unless given. PyPSA refuses one that differs by scenario (`consistency.py:1187`)
    dims: [store]
  Store_retention:
    description: share of energy kept over a snapshot — PyPSA's `(1 - standing_loss) ** elapsed hours`,
      data prep
    dims: [scenario, snapshot, store]
  Store_e_initial:
    description: energy held before the first snapshot
    dims: [scenario, store]
  Store_e_cyclic:
    description: whether the horizon closes on itself instead of opening on the initial energy
    dims: [scenario, store]
    dtype: bool
  Store_e_cyclic_per_period:
    description: whether each investment period closes on itself instead of carrying its energy on to
      the next; it overrides `e_cyclic` and `e_initial_per_period`. PyPSA reads it only under `multi_investment_periods`,
      so data prep feeds false otherwise
    dims: [scenario, store]
    dtype: bool
  Store_e_initial_per_period:
    description: whether each investment period opens on the initial energy instead of carrying the previous
      period's; PyPSA reads it only under `multi_investment_periods`, so data prep feeds false otherwise
    dims: [scenario, store]
    dtype: bool
  Store_opens_late:
    description: whether a snapshot is the first a store stands in, where that is not the first of the
      horizon — PyPSA's `active.cumsum() == 1` over the snapshots it stands in, past the first snapshot,
      data prep; false in a run where every store stands throughout
    dims: [snapshot, store]
    dtype: bool
  Store_inactive_snapshots:
    description: how many snapshots a store does not stand in — PyPSA's `(~active).sum()`, data prep.
      A cyclic store reaches back this many snapshots further, so it closes on the last snapshot it stands
      in
    dims: [store]
    dtype: int
  Store_marginal_cost:
    description: cost of one unit of power delivered
    dims: [scenario, snapshot, store]
  Store_marginal_cost_quadratic:
    description: cost of the square of the net power delivered, so charging costs as much as delivering
    dims: [scenario, snapshot, store]
  Store_marginal_cost_storage:
    description: cost of one unit of energy held over one snapshot
    dims: [scenario, snapshot, store]
variables:
  Generator_p:
    description: '`Generator-p` — output of a generator in a snapshot'
    dims: [scenario, snapshot, generator]
    where: Generator_active
  Store_e:
    description: '`Store-e` — energy held at the end of a snapshot'
    dims: [scenario, snapshot, store]
    where: Store_active
  Store_p:
    description: '`Store-p` — power delivered to the bus; charging is negative'
    dims: [scenario, snapshot, store]
    where: Store_active
  Store_e_nom_ext:
    description: '`Store-e_nom` — nominal capacity where it is a decision; the parameter of the same PyPSA
      name carries the fixed regime'
    dims: [store]
    where: Store_e_nom_extendable
  CVaR_a:
    description: '`CVaR-a` — how far a scenario''s operating cost exceeds the tail''s start; nothing where
      it does not'
    dims: [scenario]
    bounds: {lower: 0}
  CVaR_theta:
    description: '`CVaR-theta` — where the tail starts, the value at risk'
    dims: []
  CVaR:
    description: '`CVaR` — the tail''s average cost, what the objective prices at `omega`'
    dims: []
constraints:
  Generator_fix_p_lower:
    description: '`Generator-fix-p-lower` — a fixed generator outputs at least its minimum'
    dims: [scenario, snapshot, generator]
    where: not Generator_p_nom_extendable AND not Generator_committable AND Generator_active
    expression: Generator_p >= Generator_p_min_pu * Generator_p_nom
  Generator_fix_p_upper:
    description: '`Generator-fix-p-upper` — a fixed generator outputs at most what is available'
    dims: [scenario, snapshot, generator]
    where: not Generator_p_nom_extendable AND not Generator_committable AND Generator_active
    expression: Generator_p <= Generator_p_max_pu * Generator_p_nom
  Store_ext_e_lower:
    description: '`Store-ext-e-lower` — an extendable store holds at least its floor of the chosen build'
    dims: [scenario, snapshot, store]
    where: Store_e_nom_extendable AND Store_active
    expression: Store_e >= Store_e_min_pu * Store_e_nom_ext
  Store_ext_e_upper:
    description: '`Store-ext-e-upper` — an extendable store holds at most the chosen build'
    dims: [scenario, snapshot, store]
    where: Store_e_nom_extendable AND Store_active
    expression: Store_e <= Store_e_max_pu * Store_e_nom_ext
  Store_ext_e_nom_lower:
    description: '`Store-ext-e_nom-lower` — the chosen build is at least its floor in every scenario'
    dims: [scenario, store]
    where: Store_e_nom_extendable
    expression: Store_e_nom_ext >= Store_e_nom_min
  Store_ext_e_nom_upper:
    description: '`Store-ext-e_nom-upper` — the chosen build is at most its cap in every scenario; a cap
      of infinity is no row'
    dims: [scenario, store]
    where: Store_e_nom_extendable AND Store_e_nom_max
    expression: Store_e_nom_ext <= Store_e_nom_max
  Store_energy_balance:
    description: '`Store-energy_balance` — the energy carried in, less what is delivered to the bus'
    dims: [scenario, snapshot, store]
    where: Store_active
    expression: Store_e == Store_energy_carried_in - Store_p * snapshot_weightings_stores
  Bus_nodal_balance:
    description: '`Bus-nodal_balance` — what is generated at a bus, storage dispatch and stores included,
      less what the links take away, plus what arrives over them after losses and any delay at every port
      they deliver to, each process port drawing or delivering at its own rate and each passive branch
      carrying its flow, meets the load there, less half of every incident line''s and transformer''s
      loss — PyPSA dissipates a branch''s loss half at either end. Each generator, storage unit, store
      and load term enters with its component''s `sign` (`constraints.py:1428-1429`, `:1538`), and an
      inactive load not at all. A bus nothing is attached to has no row; PyPSA refuses one that carries
      load, and this file does not yet.'
    dims: [scenario, snapshot, bus]
    expression: Bus_injection == 0
  Carrier_growth_limit:
    description: '`Carrier-growth_limit` — what a carrier adds across its extendable components in a period,
      counting each build in the first period it stands in, is at most its allowance plus a share of what
      it added the period before; the first period has no predecessor, so `edge=0` leaves it the bare
      allowance'
    dims: [carrier, period]
    where: Carrier_max_growth
    expression: Carrier_additions - shift(Carrier_additions, along=period, offset=1, edge=0) * Carrier_relative_growth
      <= Carrier_max_growth
expressions:
  Store_energy_carried_in:
    description: the energy a store opens a snapshot with — at the first snapshot it stands in, its last
      such snapshot's less standing loss where it is cyclic and the given initial energy, which no standing
      loss has touched yet, where it is not; the previous snapshot's less standing loss otherwise. A store
      built in a later period opens in that period, and a cyclic one that retires closes on its own last
      snapshot. Per period, the same holds with each investment period as the horizon
    dims: [scenario, snapshot, store]
    cases:
      cyclic: {when: Store_e_cyclic AND NOT Store_e_cyclic_per_period AND NOT Store_e_initial_per_period
          AND (position(snapshot) == 0 OR Store_opens_late), expression: 'Store_retention * shift(shift(Store_e,
          along=snapshot, offset=1, edge=''wrap''), along=snapshot, offset=Store_inactive_snapshots, edge=''wrap'')'}
      opening: {when: NOT Store_e_cyclic AND NOT Store_e_cyclic_per_period AND NOT Store_e_initial_per_period
          AND (position(snapshot) == 0 OR Store_opens_late), expression: Store_e_initial}
      period_cyclic: {when: Store_e_cyclic_per_period, expression: 'Store_retention * shift(Store_e, along=snapshot,
          offset=1, edge=''wrap'', by=snapshot_period, within=period)'}
      period_opening: {when: 'Store_e_initial_per_period AND NOT Store_e_cyclic_per_period AND position(snapshot,
          by=snapshot_period, within=period) == 0', expression: Store_e_initial}
    otherwise: Store_retention * shift(Store_e, along=snapshot, offset=1)
  total_cost:
    dims: []
    expression: Store_capex + risk_weighted_opex
    description: what the system costs — capacity once per active period at its expected cost over the
      scenarios, operation in expectation over the scenarios, and a share of it at the tail
  Carrier_additions:
    dims: [period, carrier]
    expression: Store_additions
    description: what a carrier adds in a period — every extendable component of that carrier, counting
      each build in the first period it stands in. Like PyPSA, it sums only the components that carry
      a carrier attribute, so a transformer, which has none, counts in no carrier
  Carrier_relative_growth:
    description: the share of the previous period's additions a carrier's growth limit reads — PyPSA's
      `max_relative_growth` clipped at zero, so a negative share adds nothing and never tightens the limit
    dims: [carrier]
    cases:
      positive: {when: Carrier_max_relative_growth > 0, expression: Carrier_max_relative_growth}
    otherwise: 0
  Bus_injection:
    dims: [scenario, snapshot, bus]
    expression: (Generator_injection + Load_injection) + Store_injection
    description: what every component puts into a bus, less what it takes out of it; PyPSA writes each
      term into the balance, and a load on its right-hand side
  Store_capex: {expression: sum(scenario_weight * Store_e_nom_ext * Store_capital_cost * Store_capital_weight)}
  risk_weighted_opex: {expression: '(1 - CVaR_omega) * sum(scenario_weight * scenario_opex, over=scenario)
      + CVaR_omega * CVaR'}
  Store_additions: {expression: 'sum(Store_e_nom_ext * Store_first_active, by=Store_carrier, over=store,
      into=carrier)'}
  Generator_injection: {expression: 'sum(Generator_sign * Generator_p, by=Generator_bus, over=generator,
      into=bus)'}
  Load_injection: {expression: 'sum(Load_demand, by=Load_bus, over=load, into=bus)'}
  Store_injection: {expression: 'sum(Store_sign * Store_p, by=Store_bus, over=store, into=bus)'}
  scenario_opex:
    dims: [scenario]
    expression: Generator_opex + Store_opex
    description: what a future costs to run — every operating term, weighted by the snapshot's hours and
      its period, before the scenario's own weight; a start and a stop cost what they cost, unweighted,
      as PyPSA adds them (`optimize.py:414-429`)
  Load_demand:
    description: what a load draws from its bus's balance — its demand times its sign where it is active,
      nothing where it is not, since PyPSA drops an inactive load from the balance (`constraints.py:1537-1538`)
    dims: [scenario, snapshot, load]
    cases:
      active: {when: Load_active, expression: Load_sign * Load_p_set}
    otherwise: 0
  Generator_opex: {expression: 'sum(sum(((Generator_p * Generator_marginal_cost) * snapshot_weightings_objective)
      * at(period_weight_objective, by=snapshot_period, over=period, into=snapshot), over=generator),
      over=snapshot) + sum(sum((((Generator_p * Generator_p) * Generator_marginal_cost_quadratic) * snapshot_weightings_objective)
      * at(period_weight_objective, by=snapshot_period, over=period, into=snapshot), over=generator),
      over=snapshot)'}
  Store_opex: {expression: 'sum(sum(((Store_p * Store_marginal_cost) * snapshot_weightings_objective)
      * at(period_weight_objective, by=snapshot_period, over=period, into=snapshot), over=store), over=snapshot)
      + sum(sum((((Store_p * Store_p) * Store_marginal_cost_quadratic) * snapshot_weightings_objective)
      * at(period_weight_objective, by=snapshot_period, over=period, into=snapshot), over=store), over=snapshot)
      + sum(sum(((Store_e * Store_marginal_cost_storage) * snapshot_weightings_objective) * at(period_weight_objective,
      by=snapshot_period, over=period, into=snapshot), over=store), over=snapshot)'}
objective: {sense: minimize, expression: total_cost}

from differential.pypsa.prep import relation, static, varying, weighting


n = build()  # the network from the PyPSA tab

sources = {
    'snapshot': pl.Series('snapshot', list(timesteps(n)), dtype=pl.Datetime('us')),
    'bus': pl.Series('bus', list(names(n.buses.index).astype(str)), dtype=pl.String),
        **{
            dim: pl.Series(dim, list(names(n.static(component).index).astype(str)), dtype=pl.String)
            for component, dim in DIM.items()
        },
        **scenarios(n),
        **periods(n),
        **carriers(n, multi),
    'Generator_bus': relation(n, 'Generator', 'bus'),
    'Load_bus': relation(n, 'Load', 'bus'),
    'Store_bus': relation(n, 'Store', 'bus'),
    'snapshot_weightings_objective': weighting(n, 'objective'),
    'Generator_sign': per_component('Generator', first_scenario(n.generators['sign'])),
    'Load_p_set': varying(n, 'Load', 'p_set'),
    'Load_sign': per_component('Load', first_scenario(loads['sign'])),
    'Load_active': per_component('Load', first_scenario(loads['active']), bool),
    'snapshot_weightings_stores': weighting(n, 'stores'),
}

with sps.solve('differential/pypsa/rungs/rung_39_negative_relative_growth.yaml', sources) as solution:
    solution.objective  # 8600.0

# SPDX-FileCopyrightText: mathspec Contributors
#
# SPDX-License-Identifier: MIT

"""Rung 39: a negative relative growth adds nothing to a carrier's growth limit — PyPSA clips it at zero."""

from __future__ import annotations

from datetime import datetime

import pandas as pd

OPTIMIZE = {'multi_investment_periods': True}


def build():
    """A whole network, not the spine: a battery carrier with `max_relative_growth=-0.5` builds in both periods."""
    import pypsa

    n = pypsa.Network()
    n.snapshots = pd.MultiIndex.from_tuples(
        [(2020, datetime(2020, 1, 1, t)) for t in range(2)] + [(2030, datetime(2030, 1, 1, t)) for t in range(2)]
    )
    n.investment_periods = [2020, 2030]
    n.investment_period_weightings['objective'] = [1.0, 0.5]
    n.investment_period_weightings['years'] = [10.0, 10.0]
    n.snapshot_weightings['objective'] = [2.0, 1.5, 2.5, 2.0]
    n.add('Bus', 'grid')
    n.add('Carrier', 'solar')
    n.add('Carrier', 'gas')
    n.add('Carrier', 'battery', max_growth=20, max_relative_growth=-0.5)
    n.add('Generator', 'solar', bus='grid', carrier='solar', p_nom=100, marginal_cost=1, p_max_pu=[1, 0, 1, 0])
    n.add('Generator', 'backup', bus='grid', carrier='gas', p_nom=200, marginal_cost=80)
    n.add(
        'Store',
        'tank20',
        bus='grid',
        carrier='battery',
        e_nom_extendable=True,
        e_nom_max=100,
        capital_cost=10,
        build_year=2020,
        lifetime=30,
    )
    n.add(
        'Store',
        'tank30',
        bus='grid',
        carrier='battery',
        e_nom_extendable=True,
        e_nom_max=100,
        capital_cost=8,
        build_year=2030,
        lifetime=30,
    )
    n.add('Load', 'town', bus='grid', p_set=[40, 60, 40, 80])
    return n

n = build()
n.optimize(solver_name='highs')
n.objective  # 8600.0