Rung 10: quadratic costs — a marginal cost quadratic in output, stated by pypsa_quadratic.yaml¶

description: 'The quadratic class of a plain `n.optimize()`: PyPSA''s `marginal_cost_quadratic`, stated
  on rung 1''s transport surface in a file of its own. One file cannot carry a quadratic objective beside
  commitment''s integer variables and still solve on HiGHS, because degree is the spec''s property and
  not the data''s. So the class a free solver takes as a QP lives here, and `examples/pypsa.yaml` stays
  the mixed-integer one. PyPSA also carries the attribute on storage units and stores; each is one more
  term of the same shape.'
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'}
  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:
  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]
  Generator_marginal_cost_quadratic:
    description: cost of the square 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]
  Link_marginal_cost_quadratic:
    description: cost of the square 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]
  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]
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
  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)
objective: {sense: minimize, description: 'operating cost, linear and quadratic, each snapshot weighted
    by the hours it stands for', expression: sum(Generator_p * Generator_marginal_cost * snapshot_weightings_objective)
    + sum(Generator_p * Generator_p * Generator_marginal_cost_quadratic * snapshot_weightings_objective)
    + sum(Link_p * Link_marginal_cost * snapshot_weightings_objective) + sum(Link_p * Link_p * Link_marginal_cost_quadratic
    * snapshot_weightings_objective)}

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_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'),
    'Generator_marginal_cost_quadratic': varying(n, 'Generator', 'marginal_cost_quadratic'),
    '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'),
    'Link_marginal_cost_quadratic': varying(n, 'Link', 'marginal_cost_quadratic'),
    'Load_p_set': varying(n, 'Load', 'p_set'),
}

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

"""Rung 10: quadratic costs — a marginal cost quadratic in output, stated by `pypsa_quadratic.yaml`."""

from __future__ import annotations

import spine

#: This rung binds a file of its own.
MODEL = 'pypsa_quadratic.yaml'


def build():
    """The spine plus this rung's additions, as a ``pypsa.Network``."""
    n = spine.build()
    n.add('Bus', 'village')
    n.add('Generator', 'steam', bus='north', p_nom=80, marginal_cost=5, marginal_cost_quadratic=0.08)
    n.add('Generator', 'engine', bus='north', p_nom=80, marginal_cost=20, marginal_cost_quadratic=0.01)
    n.add(
        'Link',
        'wire2',
        bus0='north',
        bus1='village',
        p_nom=40,
        p_min_pu=-1,
        efficiency=0.9,
        marginal_cost=1,
        marginal_cost_quadratic=0.02,
    )
    n.add('Load', 'village_load', bus='village', p_set=15)
    n.add('Load', 'extra10', bus='north', p_set=[30, 50, 40, 60])
    return n

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