PyPSA LOPF — two coordinates on one dimension¶

description: >-
  PyPSA linear optimal power flow on a meshed AC-DC network whose
  generators sit on a bus and burn a carrier, with capacity to build and a CO2
  budget priced through the second map: the nodal balance groups through the
  bus coordinate, and the budget reads an emissions rate back down through the
  carrier. PyPSA's own ac-dc-meshed example.
  Optimum 18441021.477729216, from PyPSA itself.

dimensions:
  snapshot:
    description: dispatch periods
    dtype: int
  bus:
    description: network nodes
    dtype: str
  carrier:
    description: what a generator burns, and what its emissions are a property of
    dtype: str
  generator:
    description: >-
      generating units, each sitting on a bus and burning a carrier — two
      coordinates on one dimension, landing on two different axes
    dtype: str
  line:
    description: passive AC lines, each joining two buses
    dtype: str
  link:
    description: controllable connections, each joining two buses
    dtype: str
  cycle:
    description: one independent loop per meshed sub-network
    dtype: str

relations:
  gen_bus:
    description: the bus a generator sits on
    key: generator
    values: bus
  gen_carrier:
    description: the carrier a generator burns
    key: generator
    values: carrier
  line_from:
    description: the bus a line leaves
    key: line
    values: bus
  line_to:
    description: the bus a line arrives at
    key: line
    values: bus
  link_from:
    description: the bus a link leaves
    key: link
    values: bus
  link_to:
    description: the bus a link arrives at
    key: link
    values: bus

parameters:
  load:
    description: demand at each bus in each snapshot
    dims: [snapshot, bus]
  p_max_pu:
    description: share of built capacity a generator can produce in a snapshot
    dims: [snapshot, generator]
  p_nom_min:
    description: capacity a generator already has, and cannot fall below
    dims: [generator]
  marginal_cost:
    description: cost of one unit of output
    dims: [generator]
  gen_capital_cost:
    description: annualised cost of a unit of generator capacity
    dims: [generator]
  efficiency:
    description: share of the carrier's energy a generator turns into output
    dims: [generator]
  co2_per_mwh:
    description: emissions per unit of carrier burned, a property of the carrier
    dims: [carrier]
  line_capital_cost:
    description: annualised cost of a unit of line capacity
    dims: [line]
  link_capital_cost:
    description: annualised cost of a unit of link capacity
    dims: [link]
  link_p_max_pu:
    description: share of its capacity a link may carry forwards
    dims: [link]
  link_p_min_pu:
    description: share of its capacity a link may carry backwards, negative by convention
    dims: [link]
  cycle_incidence:
    description: >-
      the cycle basis, as a sparse table of impedance times direction. A line
      may belong to several cycles, so this cannot be a coordinate.
    dims: [cycle, line]
  co2_limit:
    description: emissions the whole horizon is allowed
    dims: []

variables:
  p:
    description: output of a generator in a snapshot
    dims: [snapshot, generator]
    bounds:
      lower: 0
  p_nom:
    description: generator capacity to hold, built on top of what already stands
    dims: [generator]
    bounds:
      lower: p_nom_min
  f:
    description: >-
      flow on a line, signed towards its `line_to` bus — not chosen, but whatever
      the voltage law leaves
    dims: [snapshot, line]
  s_nom:
    description: line capacity to build
    dims: [line]
    bounds:
      lower: 0
  g:
    description: >-
      flow on a link, signed towards the bus it delivers at — chosen, which is
      what makes it a link and not a line
    dims: [snapshot, link]
  link_p_nom:
    description: link capacity to build
    dims: [link]
    bounds:
      lower: 0

constraints:
  within_capacity:
    description: a generator produces no more than the built capacity available to it
    dims: [snapshot, generator]
    expression: p <= p_nom * p_max_pu

  line_upper:
    dims: [snapshot, line]
    expression: f <= s_nom
  line_lower:
    dims: [snapshot, line]
    expression: f >= -s_nom
  link_upper:
    dims: [snapshot, link]
    expression: g <= link_p_nom * link_p_max_pu
  link_lower:
    dims: [snapshot, link]
    expression: g >= link_p_nom * link_p_min_pu

  nodal_balance:
    description: >-
      what is generated at a bus plus what arrives over the lines and links
      meets the load there
    dims: [snapshot, bus]
    expression: >-
      sum(p, by=gen_bus, over=generator, into=bus)
      + sum(f, by=line_to, over=line, into=bus) - sum(f, by=line_from, over=line, into=bus)
      + sum(g, by=link_to, over=link, into=bus) - sum(g, by=link_from, over=link, into=bus)
      == load

  kirchhoff_voltage_law:
    description: around each independent cycle the impedance-weighted flows sum to zero
    dims: [snapshot, cycle]
    expression: sum(f * cycle_incidence, over=line) == 0

  co2_budget:
    description: >-
      PyPSA's primary-energy constraint — a generator's emissions are its
      output divided by its efficiency, priced at its carrier's rate, and the
      horizon's total stays inside the budget
    dims: []
    expression: >-
      sum(sum(p * at(co2_per_mwh, by=gen_carrier, over=carrier, into=generator) / efficiency, over=generator), over=snapshot)
      <= co2_limit

objective:
  sense: minimize
  description: what the fleet costs to run, plus what the generation and network capacity cost to build
  expression: >-
    sum(p * marginal_cost)
    + sum(p_nom * gen_capital_cost)
    + sum(s_nom * line_capital_cost)
    + sum(link_p_nom * link_capital_cost)

# sources: parameter name -> frame or parquet path
with sps.solve('examples/ports/pypsa_ac_dc.yaml', sources) as solution:
    solution.objective  # 18441021.477729216
    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``.

    The port carries its cycle basis as ``cycle_incidence`` because computing
    one is a graph algorithm and so data preparation; PyPSA derives its own
    from ``line_x`` / ``line_r``, which is why both are in the instance. The
    two must describe the same cycle space, and the objectives agreeing is
    what says they do.
    """
    n = pypsa.Network()
    n.set_snapshots(tables['snapshot']['snapshot'].tolist())

    carrier = _series(tables['bus_carrier'], 'bus')
    for bus, ct in carrier.items():
        n.add('Bus', bus, carrier=ct)

    co2 = _series(tables['co2_per_mwh'], 'carrier')
    for name, rate in co2.items():
        n.add('Carrier', name, co2_emissions=rate)

    generators = tables['generator'].set_index('generator')
    p_max_pu = _wide(tables['p_max_pu'], 'generator')
    for g, row in generators.iterrows():
        n.add(
            'Generator',
            g,
            bus=row['bus'],
            carrier=row['carrier'],
            p_nom_extendable=True,
            p_nom=_series(tables['gen_p_nom_existing'], 'generator')[g],
            p_nom_min=_series(tables['p_nom_min'], 'generator')[g],
            marginal_cost=_series(tables['marginal_cost'], 'generator')[g],
            capital_cost=_series(tables['gen_capital_cost'], 'generator')[g],
            efficiency=_series(tables['efficiency'], 'generator')[g],
            p_max_pu=p_max_pu[g],
        )

    for line, ends in tables['line'].set_index('line').iterrows():
        bus0, bus1 = ends['line_from'], ends['line_to']
        n.add(
            'Line',
            line,
            bus0=bus0,
            bus1=bus1,
            x=_series(tables['line_x'], 'line')[line],
            r=_series(tables['line_r'], 'line')[line],
            s_nom=_series(tables['line_s_nom_existing'], 'line')[line],
            s_nom_extendable=True,
            capital_cost=_series(tables['line_capital_cost'], 'line')[line],
        )

    for link, ends in tables['link'].set_index('link').iterrows():
        bus0, bus1 = ends['link_from'], ends['link_to']
        n.add(
            'Link',
            link,
            bus0=bus0,
            bus1=bus1,
            p_nom_extendable=True,
            p_nom=_series(tables['link_p_nom_existing'], 'link')[link],
            p_min_pu=_series(tables['link_p_min_pu'], 'link')[link],
            p_max_pu=_series(tables['link_p_max_pu'], 'link')[link],
            capital_cost=_series(tables['link_capital_cost'], 'link')[link],
        )

    load = _wide(tables['load'], 'bus')
    for bus in load.columns:
        if load[bus].any():
            n.add('Load', bus, bus=bus, p_set=load[bus])

    n.add(
        'GlobalConstraint',
        'co2_limit',
        type='primary_energy',
        carrier_attribute='co2_emissions',
        sense='<=',
        constant=float(tables['co2_limit']),
    )
    return n