Rung 15: two investment periods — build years, lifetimes, period weights and a carrier's growth limit, stated by pypsa_multi_period.yaml¶

description: 'The multi-period class of a plain `n.optimize()`: `multi_investment_periods`, stated on
  rungs 1 and 3 in a file of its own. A snapshot belongs to an investment period, an asset stands in the
  periods its build year and lifetime span, and capacity is paid once per period it stands in, each period
  weighted; a carrier may grow only so much per period. Which snapshots an asset is active in is data
  prep, because a `where` reaches only the frame''s own dimensions. A dimension a run may not have cannot
  ride on `examples/pypsa.yaml`, so this class lives here.'
dimensions:
  snapshot: {description: 'dispatch periods, positions across every investment period', dtype: datetime}
  period: {description: investment periods — PyPSA's `investment_periods`, dtype: int}
  bus: {description: network nodes}
  generator: {description: 'generating units, each on one bus'}
  link: {description: 'controllable connections, each from one bus to the buses it delivers to'}
  link_output: {description: 'a link''s output ports, one label per port a link declares — PyPSA''s `bus1`,
      `bus2`, … columns read long, so a link of any number of output ports is one term in the balance,
      data prep'}
  load: {description: 'demands, each on one bus'}
  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_carrier: {description: the carrier a generator converts from, key: generator, values: carrier}
  Generator_bus: {description: the bus a generator sits on, key: generator, values: bus}
  Link_bus0: {description: the bus a link leaves, key: link, values: bus}
  Link_output_link: {description: the link an output port belongs to, key: link_output, values: link}
  Link_output_bus: {description: 'the bus an output port delivers to — PyPSA''s `bus1`, `bus2`, … columns.
      A link of three output ports is three labels here rather than a third relation, so the file states
      any number of them', key: link_output, values: bus}
  Load_bus: {description: the bus a load sits on, key: load, values: bus}
parameters:
  snapshot_weightings_objective:
    description: PyPSA's `snapshot_weightings.objective` — hours a snapshot stands for in the cost
    dims: [snapshot]
  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
  Generator_capital_weight:
    description: the sum of period weights a generator stands in — PyPSA's `active * period_weighting`,
      summed, data prep
    dims: [generator]
  Generator_first_active:
    description: one in the first period a generator stands in, zero elsewhere — PyPSA's `active.cumsum()
      == 1`, data prep
    dims: [period, generator]
  Carrier_max_growth:
    description: most capacity of a carrier that may be added in a period; no value means no limit
    dims: [carrier]
  Carrier_max_relative_growth:
    description: share of the previous period's additions that may be added on top
    dims: [carrier]
  Generator_p_nom:
    description: nominal power
    dims: [generator]
  Generator_p_nom_extendable:
    description: whether the nominal power is a decision
    dims: [generator]
    dtype: bool
  Generator_p_nom_min:
    description: least nominal power an extendable generator may be built at
    dims: [generator]
  Generator_p_nom_max:
    description: most nominal power an extendable generator may be built at
    dims: [generator]
  Generator_capital_cost:
    description: cost of one unit of nominal power — PyPSA's `capital_cost`, periodized as an annuity
      in data prep
    dims: [generator]
  Generator_p_min_pu:
    description: least output, per unit of nominal power
    dims: [snapshot, generator]
  Generator_p_max_pu:
    description: most output, per unit of nominal power — an availability profile
    dims: [snapshot, generator]
  Generator_marginal_cost:
    description: cost of one unit of output
    dims: [snapshot, generator]
  Link_p_nom:
    description: nominal power
    dims: [link]
  Link_p_min_pu:
    description: least flow, per unit of nominal power — negative for a link that carries both ways
    dims: [snapshot, link]
  Link_p_max_pu:
    description: most flow, per unit of nominal power
    dims: [snapshot, link]
  Link_efficiency:
    description: share of the flow that arrives at an output port, PyPSA's `efficiency`, `efficiency2`,
      … read long — negative where that port consumes rather than delivers
    dims: [link_output]
  Link_marginal_cost:
    description: cost of one unit of flow
    dims: [snapshot, link]
  Load_p_set:
    description: demand
    dims: [snapshot, load]
variables:
  Generator_p:
    description: '`Generator-p` — output of a generator in a snapshot'
    dims: [snapshot, generator]
    where: Generator_active
  Link_p:
    description: '`Link-p` — PyPSA''s `p0`, the flow measured at the `Link_bus0` end: a positive value
      withdraws there and injects at every bus the link''s output ports deliver to'
    dims: [snapshot, link]
  Generator_p_nom_ext:
    description: '`Generator-p_nom` — nominal power where it is a decision; the parameter of the same
      PyPSA name carries the fixed regime'
    dims: [generator]
    where: Generator_p_nom_extendable
constraints:
  Generator_fix_p_lower:
    description: '`Generator-fix-p-lower` — a generator outputs at least its minimum'
    dims: [snapshot, generator]
    where: not Generator_p_nom_extendable AND Generator_active
    expression: Generator_p >= Generator_p_min_pu * Generator_p_nom
  Generator_fix_p_upper:
    description: '`Generator-fix-p-upper` — a generator outputs at most what is available'
    dims: [snapshot, generator]
    where: not Generator_p_nom_extendable AND Generator_active
    expression: Generator_p <= Generator_p_max_pu * Generator_p_nom
  Generator_ext_p_lower:
    description: '`Generator-ext-p-lower` — an extendable generator outputs at least its minimum of the
      chosen build'
    dims: [snapshot, generator]
    where: Generator_p_nom_extendable AND Generator_active
    expression: Generator_p >= Generator_p_min_pu * Generator_p_nom_ext
  Generator_ext_p_upper:
    description: '`Generator-ext-p-upper` — an extendable generator outputs at most what is available
      of the chosen build'
    dims: [snapshot, generator]
    where: Generator_p_nom_extendable AND Generator_active
    expression: Generator_p <= Generator_p_max_pu * Generator_p_nom_ext
  Generator_ext_p_nom_lower:
    description: '`Generator-ext-p_nom-lower` — the chosen build is at least its floor'
    dims: [generator]
    where: Generator_p_nom_extendable
    expression: Generator_p_nom_ext >= Generator_p_nom_min
  Generator_ext_p_nom_upper:
    description: '`Generator-ext-p_nom-upper` — the chosen build is at most its cap; a cap of infinity
      is no row'
    dims: [generator]
    where: Generator_p_nom_extendable AND Generator_p_nom_max
    expression: Generator_p_nom_ext <= Generator_p_nom_max
  Link_fix_p_lower:
    description: '`Link-fix-p-lower` — a link carries at least its minimum, negative for the other way'
    dims: [snapshot, link]
    expression: Link_p >= Link_p_min_pu * Link_p_nom
  Link_fix_p_upper:
    description: '`Link-fix-p-upper` — a link carries at most its nominal power'
    dims: [snapshot, link]
    expression: Link_p <= Link_p_max_pu * Link_p_nom
  Bus_nodal_balance:
    description: '`Bus-nodal_balance` — what is generated at a bus, less what the links take away, plus
      what arrives over them after losses, meets the load there'
    dims: [snapshot, bus]
    expression: sum(Generator_p, by=Generator_bus, over=generator, into=bus) - sum(Link_p, by=Link_bus0,
      over=link, into=bus) + sum(at(Link_p, by=Link_output_link, over=link, into=link_output) * Link_efficiency,
      by=Link_output_bus, over=link_output, into=bus) == sum(Load_p_set, by=Load_bus, over=load, into=bus)
  Carrier_growth_limit:
    description: '`Carrier-growth_limit` — what a carrier adds 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: sum(Generator_p_nom_ext * Generator_first_active, by=Generator_carrier, over=generator,
      into=carrier) - shift(sum(Generator_p_nom_ext * Generator_first_active, by=Generator_carrier, over=generator,
      into=carrier), along=period, offset=1, edge=0) * Carrier_max_relative_growth <= Carrier_max_growth
objective: {sense: minimize, description: 'operating cost by weighted snapshot and weighted period, and
    capacity once per period it stands in', expression: 'sum(Generator_p * Generator_marginal_cost * snapshot_weightings_objective
    * at(period_weight_objective, by=snapshot_period, over=period, into=snapshot)) + sum(Link_p * Link_marginal_cost
    * snapshot_weightings_objective * at(period_weight_objective, by=snapshot_period, over=period, into=snapshot))
    + sum(Generator_p_nom_ext * Generator_capital_cost * Generator_capital_weight)'}

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


def _link_ports(n: pypsa.Network) -> pd.DataFrame:
    """A link's output ports read long — one row per port a link declares, carrying the link, the bus it delivers to and its efficiency.

    PyPSA spells the ports across columns — ``bus1``/``efficiency``, ``bus2``/``efficiency2``, … — and a
    link declares a port by naming a bus in one, so a link of any port count is as many rows here and
    one term in the balance. The label is the link and the column the port came from.
    """
    links = n.static('Link')
    blank = pd.Series('', index=links.index, dtype=str)
    frames = []
    for port in ['1', *n.components.links.additional_ports]:
        suffix = '' if port == '1' else port
        buses = links.get(f'bus{port}', blank).astype(str)
        # `efficiency`, `delay` and `cyclic_delay` are PyPSA's unsuffixed attributes: port 1
        # spells them bare and every port after it takes the number
        efficiencies = links.get(f'efficiency{suffix}', pd.Series(1.0, index=links.index)).astype(float)
        delays = links.get(f'delay{suffix}', pd.Series(0, index=links.index)).fillna(0).astype(int)
        cyclic = links.get(f'cyclic_delay{suffix}', pd.Series(False, index=links.index)).fillna(False).astype(bool)
        frame = pd.DataFrame(
            keyed(links.index, 'link')
            | {
                'bus': buses.to_numpy(),
                'value': efficiencies.to_numpy(),
                'delay': delays.to_numpy(),
                'cyclic_delay': cyclic.to_numpy(),
                'port': int(port),
            }
        )
        frames.append(frame[buses.to_numpy() != ''])
    ports = pd.concat(frames, ignore_index=True).sort_values(['link', 'port'], kind='stable')
    ports['link_output'] = ports['link'] + '_bus' + ports['port'].astype(str)
    return ports.drop(columns='port').reset_index(drop=True)


def _per_port(n: pypsa.Network, column: str, as_name: str | None = None) -> pd.DataFrame:
    """One column of the long port table keyed by ``link_output`` — what a port names, or what it carries.

    *as_name* is what the file calls it: a relation keeps its target dimension's
    own name, and every parameter over the ports lands under ``value``.
    """
    ports = _link_ports(n)
    keys = [key for key in ('scenario', 'link_output') if key in ports.columns]
    return ports[[*keys, column]].rename(columns={column: as_name or column})


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),
    'generator': pl.Series('generator', list(names(generators.index).astype(str)), dtype=pl.String),
    'link': pl.Series('link', list(names(links.index).astype(str)), dtype=pl.String),
    'link_output': pl.Series('link_output', list(pd.unique(_link_ports(n)['link_output'])), dtype=pl.String),
    'load': pl.Series('load', list(names(loads.index).astype(str)), dtype=pl.String),
    'Generator_bus': relation(n, 'Generator', 'bus'),
    'Link_bus0': relation(n, 'Link', 'bus0'),
    'Link_output_link': _per_port(n, 'link'),
    'Link_output_bus': _per_port(n, 'bus'),
    'Load_bus': relation(n, 'Load', 'bus'),
    'snapshot_weightings_objective': weighting(n, 'objective'),
    'Generator_p_nom': static(n, 'Generator', 'p_nom'),
    'Generator_p_nom_extendable': static(n, 'Generator', 'p_nom_extendable'),
    'Generator_p_nom_min': static(n, 'Generator', 'p_nom_min'),
    'Generator_p_nom_max': static(n, 'Generator', 'p_nom_max'),
    'Generator_capital_cost': static(n, 'Generator', 'capital_cost'),
    'Generator_p_min_pu': varying(n, 'Generator', 'p_min_pu'),
    'Generator_p_max_pu': varying(n, 'Generator', 'p_max_pu'),
    'Generator_marginal_cost': varying(n, 'Generator', 'marginal_cost'),
    'Link_p_nom': static(n, 'Link', 'p_nom'),
    'Link_p_min_pu': varying(n, 'Link', 'p_min_pu'),
    'Link_p_max_pu': varying(n, 'Link', 'p_max_pu'),
    'Link_efficiency': _per_port(n, 'value'),
    'Link_marginal_cost': varying(n, 'Link', 'marginal_cost'),
    'Load_p_set': varying(n, 'Load', 'p_set'),
}

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

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

"""Rung 15: two investment periods — build years, lifetimes, period weights and a carrier's growth limit, stated by `pypsa_multi_period.yaml`."""

from __future__ import annotations

from datetime import datetime

import pandas as pd

MODEL = 'pypsa_multi_period.yaml'
OPTIMIZE = {'multi_investment_periods': True}


def build():
    """A whole network, not the spine: eight snapshots over two periods, a unit that retires, two wind builds capped by growth."""
    import pypsa

    n = pypsa.Network()
    n.snapshots = pd.MultiIndex.from_tuples(
        [(2020, datetime(2020, 1, 1, t)) for t in range(4)] + [(2030, datetime(2030, 1, 1, t)) for t in range(4)]
    )
    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, 2.0, 1.5, 2.5, 2.0]
    n.add('Bus', ['north', 'south'])
    n.add('Carrier', 'wind', max_growth=50, max_relative_growth=0.5)
    n.add('Carrier', 'gas')
    n.add('Generator', 'old_gas', bus='north', carrier='gas', p_nom=40, marginal_cost=30, build_year=2010, lifetime=15)
    n.add(
        'Generator',
        'wind20',
        bus='north',
        carrier='wind',
        p_nom_extendable=True,
        p_nom_max=200,
        marginal_cost=1,
        capital_cost=100,
        build_year=2020,
        lifetime=30,
        p_max_pu=[0.8, 0.6, 0.7, 0.5, 0.8, 0.6, 0.7, 0.5],
    )
    n.add(
        'Generator',
        'wind30',
        bus='south',
        carrier='wind',
        p_nom_extendable=True,
        p_nom_max=200,
        marginal_cost=1,
        capital_cost=80,
        build_year=2030,
        lifetime=30,
        p_max_pu=[0.9, 0.7, 0.6, 0.8, 0.9, 0.7, 0.6, 0.8],
    )
    n.add(
        'Generator',
        'gas30',
        bus='south',
        carrier='gas',
        p_nom_extendable=True,
        p_nom_max=200,
        marginal_cost=40,
        capital_cost=50,
        build_year=2030,
        lifetime=30,
    )
    n.add('Link', 'wire15', bus0='north', bus1='south', p_nom=60, p_min_pu=-1, efficiency=0.95)
    n.add('Load', 'town15', bus='north', p_set=[20, 30, 25, 20, 35, 45, 40, 30])
    n.add('Load', 'port15', bus='south', p_set=[10, 20, 15, 10, 30, 40, 35, 25])
    return n

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