Rung 13: transmission losses in tangent form — a loss per line, stated by pypsa_losses.yaml¶

description: 'The lossy class of a plain `n.optimize()`: `transmission_losses` in its tangent form, stated
  on rung 6''s lines in a file of its own. A line dissipates a loss its flow buys along a fan of tangents
  to the quadratic curve, half at either end — a variable and rows the keyword adds, which no `where:`
  can add to `examples/pypsa.yaml`. The fan''s slopes and offsets are data prep, one per segment.'
dimensions:
  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'}
  line: {description: 'passive branches, each between two buses, their flow set by impedance'}
  cycle: {description: 'independent cycles of the passive network graph — the cycle basis, data prep'}
  segment: {description: 'the tangents the loss curve is approximated by, PyPSA''s `segments`', dtype: int}
  load: {description: 'demands, each on one bus'}
relations:
  Generator_bus: {description: the bus a generator sits on, key: generator, values: bus}
  Line_bus0: {description: the bus a line's flow is measured at, key: line, values: bus}
  Line_bus1: {description: the bus at a line's other end, key: line, 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]
  Generator_p_nom:
    description: nominal power
    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]
  Line_s_nom:
    description: nominal apparent power
    dims: [line]
  Line_s_nom_extendable:
    description: whether the nominal apparent power is a decision
    dims: [line]
    dtype: bool
  Line_s_max_pu:
    description: most flow either way, per unit of nominal apparent power
    dims: [snapshot, line]
  Line_s_nom_min:
    description: least nominal apparent power an extendable line may be built at
    dims: [line]
  Line_s_nom_max:
    description: most nominal apparent power an extendable line may be built at
    dims: [line]
  Line_capital_cost:
    description: cost of one unit of nominal apparent power — PyPSA's `capital_cost`, periodized as an
      annuity in data prep
    dims: [line]
  Line_cycle_weight:
    description: the line's series impedance, signed by its orientation in the cycle — the cycle basis,
      data prep; a line in no cycle has no row
    dims: [line, cycle]
  Line_loss_max:
    description: the loss at a line's rating — PyPSA's `r_pu_eff * (s_max_pu * s_nom_max)**2`, data prep
    dims: [snapshot, line]
  Line_loss_slope:
    description: the slope of a tangent to the loss curve at its segment's flow — `2 * r_pu_eff * p_k`,
      data prep
    dims: [snapshot, line, segment]
  Line_loss_offset:
    description: where that tangent meets the loss axis — `loss_k - slope_k * p_k`, negative, data prep
    dims: [snapshot, line, segment]
  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]
  Line_s:
    description: '`Line-s` — PyPSA''s `p0`, the flow measured at the `Line_bus0` end: a positive value
      withdraws there and injects at `Line_bus1`'
    dims: [snapshot, line]
  Line_s_nom_ext:
    description: '`Line-s_nom` — nominal apparent power where it is a decision; the parameter of the same
      PyPSA name carries the fixed regime'
    dims: [line]
    where: Line_s_nom_extendable
  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]
  Line_loss:
    description: '`Line-loss` — what a line dissipates carrying its flow, pushed down by the cost and
      held up by the tangents'
    dims: [snapshot, line]
    bounds: {lower: 0}
constraints:
  Generator_fix_p_lower:
    description: '`Generator-fix-p-lower` — a generator outputs at least its minimum'
    dims: [snapshot, generator]
    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]
    expression: Generator_p <= Generator_p_max_pu * Generator_p_nom
  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
  Line_fix_s_lower:
    description: '`Line-fix-s-lower` — a fixed line carries at least the negative of its rating, the loss
      counted against it'
    dims: [snapshot, line]
    where: not Line_s_nom_extendable
    expression: Line_s - Line_loss >= -Line_s_max_pu * Line_s_nom
  Line_fix_s_upper:
    description: '`Line-fix-s-upper` — a fixed line carries at most its rating, loss included'
    dims: [snapshot, line]
    where: not Line_s_nom_extendable
    expression: Line_s + Line_loss <= Line_s_max_pu * Line_s_nom
  Line_ext_s_lower:
    description: '`Line-ext-s-lower` — an extendable line carries at least the negative of its rating
      of the chosen build'
    dims: [snapshot, line]
    where: Line_s_nom_extendable
    expression: Line_s - Line_loss >= -Line_s_max_pu * Line_s_nom_ext
  Line_ext_s_upper:
    description: '`Line-ext-s-upper` — an extendable line carries at most its rating of the chosen build'
    dims: [snapshot, line]
    where: Line_s_nom_extendable
    expression: Line_s + Line_loss <= Line_s_max_pu * Line_s_nom_ext
  Line_ext_s_nom_lower:
    description: '`Line-ext-s_nom-lower` — the chosen build is at least its floor'
    dims: [line]
    where: Line_s_nom_extendable
    expression: Line_s_nom_ext >= Line_s_nom_min
  Line_ext_s_nom_upper:
    description: '`Line-ext-s_nom-upper` — the chosen build is at most its cap; a cap of infinity is no
      row'
    dims: [line]
    where: Line_s_nom_extendable AND Line_s_nom_max
    expression: Line_s_nom_ext <= Line_s_nom_max
  Kirchhoff_Voltage_Law:
    description: '`Kirchhoff-Voltage-Law` — around every independent cycle the impedance-weighted flows
      sum to nothing, which is what makes the linear power flow physical rather than transport'
    dims: [snapshot, cycle]
    expression: sum(Line_s * Line_cycle_weight, over=line) == 0
  Bus_nodal_balance:
    description: '`Bus-nodal_balance` — what is generated at a bus, plus what the links and lines bring,
      meets the load there, less half of every incident line''s loss — PyPSA dissipates a branch''s loss
      half at either end'
    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(Line_s, by=Line_bus0, over=line, into=bus)
      + sum(Line_s, by=Line_bus1, over=line, into=bus) - 0.5 * sum(Line_loss, by=Line_bus0, over=line,
      into=bus) - 0.5 * sum(Line_loss, by=Line_bus1, over=line, into=bus) == sum(Load_p_set, by=Load_bus,
      over=load, into=bus)
  Line_loss_upper:
    description: '`Line-loss_upper` — a line dissipates at most the loss at its rating'
    dims: [snapshot, line]
    expression: Line_loss <= Line_loss_max
  Line_loss_tangents_forward:
    description: '`Line-loss_tangents-{k}-1` — the loss sits above every tangent to its curve for flow
      one way; PyPSA names one row per segment `k`, this block states them all over the segment dimension'
    dims: [snapshot, line, segment]
    expression: Line_loss + Line_loss_slope * Line_s >= Line_loss_offset
  Line_loss_tangents_reverse:
    description: '`Line-loss_tangents-{k}--1` — the same fan mirrored, the loss depending on the flow''s
      magnitude'
    dims: [snapshot, line, segment]
    expression: Line_loss - Line_loss_slope * Line_s >= Line_loss_offset
objective: {sense: minimize, description: 'operating cost by weighted snapshot, plus what the lines cost
    to build', expression: sum(Generator_p * Generator_marginal_cost * snapshot_weightings_objective)
    + sum(Link_p * Link_marginal_cost * snapshot_weightings_objective) + sum(Line_s_nom_ext * Line_capital_cost)}

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


def _cycle_weights(n: pypsa.Network) -> pd.DataFrame:
    """The KVL rows PyPSA itself writes — ``n.cycle_matrix(apply_weights=True)``, reactance on AC and resistance on DC, times the 1e5 PyPSA scales every cycle row by for conditioning."""
    n.determine_network_topology()
    n.calculate_dependent_values()
    cycles = n.cycle_matrix(apply_weights=True) * 1e5
    rows = [
        {'line': str(name), 'cycle': str(cycle), 'value': float(weight)}
        for (kind, name), weights in cycles.iterrows()
        for cycle, weight in weights.items()
        if kind == 'Line' and weight
    ]
    return pd.DataFrame(rows, columns=['line', 'cycle', 'value']).astype({'value': float})


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})


def _weights(gcs: pd.DataFrame, components: pd.DataFrame, dim: str, value) -> pd.DataFrame:
    """One row per (global constraint, member): *value* returns the weight, or 0/None outside the row's set."""
    rows = [
        {'global_constraint': str(label), dim: str(name), 'value': float(v)}
        for label, gc in gcs.iterrows()
        for name, component in components.iterrows()
        if (v := value(gc, component))
    ]
    return pd.DataFrame(rows, columns=['global_constraint', dim, 'value']).astype({'value': float})


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),
    'line': pl.Series('line', list(names(lines.index).astype(str)), dtype=pl.String),
    'cycle': pl.Series('cycle', list(pd.unique(tables['Line_cycle_weight']['cycle'])), dtype=pl.String),
    'load': pl.Series('load', list(names(loads.index).astype(str)), dtype=pl.String),
    'Generator_bus': relation(n, 'Generator', 'bus'),
    'Line_bus0': relation(n, 'Line', 'bus0'),
    'Line_bus1': relation(n, 'Line', 'bus1'),
    '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_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'),
    'Line_s_nom': static(n, 'Line', 's_nom'),
    'Line_s_nom_extendable': static(n, 'Line', 's_nom_extendable'),
    'Line_s_max_pu': varying(n, 'Line', 's_max_pu'),
    'Line_s_nom_min': static(n, 'Line', 's_nom_min'),
    'Line_s_nom_max': static(n, 'Line', 's_nom_max'),
    'Line_capital_cost': static(n, 'Line', 'capital_cost'),
    'Line_cycle_weight': _cycle_weights(n),
    '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_13_losses.yaml', sources) as solution:
    solution.objective  # 10805.29588

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

"""Rung 13: transmission losses in tangent form — a loss per line, stated by `pypsa_losses.yaml`."""

from __future__ import annotations

import spine

MODEL = 'pypsa_losses.yaml'
OPTIMIZE = {'transmission_losses': {'mode': 'tangents', 'segments': 2}}


def build():
    """The spine plus a 110 kV triangle of lines, one of them extendable — ohms a real line has, so the loss stays a few percent of the flow."""
    n = spine.build()
    n.add('Bus', ['a', 'b', 'c'], v_nom=110)
    n.add('Generator', 'hydro13', bus='a', p_nom=80, marginal_cost=10)
    n.add('Generator', 'diesel13', bus='b', p_nom=80, marginal_cost=50)
    n.add('Line', 'ab13', bus0='a', bus1='b', carrier='AC', x=30, r=6, s_nom=60)
    n.add('Line', 'bc13', bus0='b', bus1='c', carrier='AC', x=60, r=9.7, s_nom=60)
    n.add(
        'Line',
        'ca13',
        bus0='c',
        bus1='a',
        carrier='AC',
        x=45,
        r=6,
        s_nom=40,
        s_nom_extendable=True,
        s_nom_max=90,
        capital_cost=4,
    )
    n.add('Load', 'town13', bus='c', p_set=[35, 55, 15, 45])
    return n

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