Rung 14: two futures and a risk preference — capacity chosen once, dispatch per scenario, stated by pypsa_stochastic.yaml¶

description: 'The two-stage class of a plain `n.optimize()`: a network with scenarios, stated on rung
  1''s transport and rung 3''s expansion in a file of its own. Everything over a snapshot spans a scenario
  as well; capacity does not — it is chosen once, before the future is known — and the cost is the expectation
  over the scenarios'' weights. With a risk preference PyPSA adds the CVaR rows: an excess per scenario
  and the tail''s average, blended into the objective. A dimension a run may not have cannot ride on `examples/pypsa.yaml`,
  so this class lives here.'
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'}
  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'}
relations:
  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:
  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 the operating cost priced at the tail
      rather than in expectation
    dims: []
  CVaR_inv_tail:
    description: PyPSA's `1 / (1 - alpha)` — the tail's own probability, inverted in data prep because
      a divisor is one factor
    dims: []
  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: [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: [scenario, 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: [scenario, snapshot, load]
variables:
  Generator_p:
    description: '`Generator-p` — output of a generator in a snapshot'
    dims: [scenario, snapshot, generator]
  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: [scenario, 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
  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 generator outputs at least its minimum'
    dims: [scenario, snapshot, generator]
    where: not Generator_p_nom_extendable
    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: [scenario, snapshot, generator]
    where: not Generator_p_nom_extendable
    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: [scenario, snapshot, generator]
    where: Generator_p_nom_extendable
    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: [scenario, snapshot, generator]
    where: Generator_p_nom_extendable
    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: [scenario, 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: [scenario, 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: [scenario, 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)
  CVaR_excess:
    description: '`CVaR-excess-{s}` — a scenario''s operating cost beyond the tail''s start is its excess;
      PyPSA names one row per scenario'
    dims: [scenario]
    expression: CVaR_a - scenario_opex + CVaR_theta >= 0
  CVaR_def:
    description: '`CVaR-def` — the tail''s average is at least where it starts plus the expected excess
      over the tail''s probability'
    dims: []
    expression: CVaR_theta + CVaR_inv_tail * sum(scenario_weight * CVaR_a, over=scenario) <= CVaR
expressions:
  scenario_opex: {description: 'what a future costs to run — the operating terms, before their weight',
    expression: 'sum(sum(Generator_p * Generator_marginal_cost * snapshot_weightings_objective, over=generator),
      over=snapshot) + sum(sum(Link_p * Link_marginal_cost * snapshot_weightings_objective, over=link),
      over=snapshot)'}
objective: {sense: minimize, description: 'capacity once, operation in expectation, and a share of it
    at the tail', expression: 'sum(Generator_p_nom_ext * Generator_capital_cost) + (1 - CVaR_omega) *
    sum(scenario_weight * scenario_opex, over=scenario) + CVaR_omega * CVaR'}

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_14_stochastic.yaml', sources) as solution:
    solution.objective  # 9267.386667

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

"""Rung 14: two futures and a risk preference — capacity chosen once, dispatch per scenario, stated by `pypsa_stochastic.yaml`."""

from __future__ import annotations

import spine

MODEL = 'pypsa_stochastic.yaml'


def build():
    """The spine plus an extendable wind unit whose availability and the south's load differ between a calm and a stormy future."""
    n = spine.build()
    n.add('Generator', 'wind14', bus='south', p_nom_extendable=True, p_nom_max=100, marginal_cost=1, capital_cost=20)
    n.add('Load', 'port14', bus='south')
    n.set_scenarios({'calm': 0.6, 'stormy': 0.4})
    n.c.loads.dynamic.p_set[('calm', 'port14')] = [10, 20, 15, 10]
    n.c.loads.dynamic.p_set[('stormy', 'port14')] = [40, 60, 50, 30]
    n.c.generators.dynamic.p_max_pu[('calm', 'wind14')] = [0.9, 0.7, 0.8, 0.6]
    n.c.generators.dynamic.p_max_pu[('stormy', 'wind14')] = [0.3, 0.2, 0.4, 0.1]
    n.set_risk_preference(alpha=0.5, omega=0.3)
    return n

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