PyPSA transmission losses — a quadratic curve, held up by its own tangents¶

description: >-
  PyPSA transmission losses, tangent form: the quadratic loss on a line,
  underestimated from below by a fan of tangent half-planes and subtracted half
  at each end. No auxiliary variable and no piecewise construct — a tangent is
  one linear row per segment and per sign of the flow.
  Optimum 3692.705905599654, from PyPSA itself.

dimensions:
  snapshot:
    description: dispatch periods
    dtype: int
  bus:
    description: network nodes
    dtype: str
  generator:
    description: generating units, each sitting on one bus
    dtype: str
  line:
    description: passive branches, each joining two buses
    dtype: str
  segment:
    description: the tangent points the loss curve is approximated at
    dtype: int

relations:
  gen_bus:
    description: the bus a generator sits on
    key: generator
    values: bus
  from:
    description: the bus a line leaves
    key: line
    values: bus
  to:
    description: the bus a line arrives at
    key: line
    values: bus

parameters:
  p_nom:
    description: installed capacity of a generator
    dims: [generator]
  marginal_cost:
    description: cost of one unit of output
    dims: [generator]
  s_nom:
    description: most a line may carry towards its `to` bus
    dims: [line]
  neg_s_nom:
    description: most a line may carry the other way, negative by convention
    dims: [line]
  loss_max:
    description: >-
      the loss at a line's rating — the top of the curve being approximated,
      carried as a column because a bound takes a name or a number, and given
      only for the lines that dissipate anything
    dims: [line]
  loss_slope:
    description: >-
      the slope of this segment's half-plane — how much loss the flow buys
      along it. Where the segments come from is the instance's business, not
      the model's: a tangent to the loss curve and a secant across it both
      arrive here as a slope and an offset.
    dims: [line, segment]
  loss_offset:
    description: >-
      where this segment's half-plane meets the loss axis, negative for a curve
      through the origin
    dims: [line, segment]
  load:
    description: demand at each bus in each snapshot
    dims: [snapshot, bus]

variables:
  p:
    description: output of a generator in a snapshot
    dims: [snapshot, generator]
    bounds:
      lower: 0
      upper: p_nom
  f:
    description: >-
      flow on a line, signed towards its `to` bus — unbounded here, because the
      rating covers the flow and its loss and so is a row rather than a bound
    dims: [snapshot, line]
  loss:
    description: >-
      the energy a line dissipates carrying its flow — pushed down by the
      objective and held up by the tangents, so it settles on the approximated
      curve rather than needing an equality of its own. A line with no
      resistance dissipates nothing, which is a loss of zero rather than a
      quantity with no value, so the balances and ratings that name it keep
      their rows.
    dims: [snapshot, line]
    where: loss_max
    absence: zero
    bounds:
      lower: 0
      upper: loss_max

constraints:
  nodal_balance:
    description: >-
      what is generated at a bus plus what arrives over the lines meets the load
      there, less half of each incident line's loss — PyPSA's convention is that
      a branch dissipates half at either end
    dims: [snapshot, bus]
    expression: >-
      sum(p, by=gen_bus, over=generator, into=bus)
      + sum(f, by=to, over=line, into=bus) - sum(f, by=from, over=line, into=bus)
      - 0.5 * sum(loss, by=from, over=line, into=bus) - 0.5 * sum(loss, by=to, over=line, into=bus)
      == load

  within_rating_forward:
    description: >-
      a line's rating limits what it carries plus what it dissipates, so the
      loss eats into the capacity rather than riding on top of it
    dims: [snapshot, line]
    expression: f + loss <= s_nom

  within_rating_reverse:
    description: the same limit for flow the other way
    dims: [snapshot, line]
    expression: f - loss >= neg_s_nom

  loss_above_segment_forward:
    description: >-
      the loss sits above every one of its half-planes, which for a convex
      curve is the whole approximation. Only the lines that have a curve get a
      fan.
    dims: [snapshot, line, segment]
    where: loss_max
    expression: loss + loss_slope * f >= loss_offset

  loss_above_segment_reverse:
    description: the same fan mirrored, because the loss depends on the flow's magnitude
    dims: [snapshot, line, segment]
    where: loss_max
    expression: loss - loss_slope * f >= loss_offset

objective:
  sense: minimize
  description: what the fleet costs to run; the losses are paid for as extra generation
  expression: sum(p * marginal_cost)

# sources: parameter name -> frame or parquet path
with sps.solve('examples/ports/pypsa_losses.yaml', sources) as solution:
    solution.objective  # 24114.237385131008
    solution.dual('nodal_balance')

def build(tables: dict[str, pd.DataFrame]) -> pypsa.Network:
    """The port's tables as a PyPSA network, column for column.

    ``tables`` is the same mapping the specsolve call attaches as ``sources``.

    PyPSA is given ``r``, ``x`` and ``s_nom`` and derives the tangents itself.
    The port is given the tangents, because a slope of ``2 * r * p_k`` is
    arithmetic and the language's coefficients take a name or a number — the
    same reason ``pypsa_storage`` ships ``soc_max`` rather than a ratio. Both
    sides therefore describe one model from the same instance, and
    ``SEGMENTS`` is the one number that has to agree between them.
    """
    n = pypsa.Network()
    n.set_snapshots(tables['snapshot']['snapshot'])
    n.add('Bus', tables['bus']['bus'])

    lines: pd.DataFrame = tables['line'].set_index('line')
    n.add(
        'Line',
        lines.index,
        bus0=lines['from'],
        bus1=lines['to'],
        r=tables['r'].set_index('line')['value'],
        x=tables['x'].set_index('line')['value'],
        s_nom=tables['s_nom'].set_index('line')['value'],
    )

    generators: pd.DataFrame = tables['generator'].set_index('generator')
    n.add(
        'Generator',
        generators.index,
        bus=generators['gen_bus'],
        p_nom=tables['p_nom'].set_index('generator')['value'],
        marginal_cost=tables['marginal_cost'].set_index('generator')['value'],
    )

    load: pd.DataFrame = tables['load'].pivot(index='snapshot', columns='bus', values='value')
    for bus in load.columns:
        n.add('Load', f'load_{bus}', bus=bus, p_set=load[bus])
    return n