Rung 11: PyPSA's own ac_dc_meshed example, whole — meshed AC and DC, extendable lines, links and generators, carriers, a CO2 budget¶

description: The spec of the model a plain `n.optimize()` builds, in one file. Every declaration is named
  `Component_attribute` after the PyPSA statement it stands for, and each constraint's description opens
  with the linopy name PyPSA gives that row, so the two can be read side by side. PyPSA's regimes — extendable,
  committable — are data columns and become `where:` masks. Bounds are the explicit rows PyPSA writes,
  so their duals are row duals. Parameters no PyPSA table carries verbatim are computed in data prep and
  say so in their description.
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'}
  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'}
  global_constraint: {description: 'PyPSA''s `GlobalConstraint` rows, one label per declared limit'}
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}
  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}
parameters:
  snapshot_weightings_objective:
    description: PyPSA's `snapshot_weightings.objective` — hours a snapshot stands for in the cost
    dims: [snapshot]
  Generator_p_nom_extendable:
    description: whether the nominal power is a decision
    dims: [generator]
    dtype: bool
  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_committable:
    description: whether output is gated by an on/off status decision
    dims: [generator]
    dtype: bool
  Link_p_nom_extendable:
    description: whether the nominal power is a decision
    dims: [link]
    dtype: bool
  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_output_delay:
    description: snapshots a port's delivery lags its link's flow — PyPSA's `delay`, `delay2`, … read
      long, in `snapshot_weightings.generators` units, which the file states as whole snapshots; zero
      for a port that delivers at once
    dims: [link_output]
    dtype: int
  Link_output_cyclic_delay:
    description: whether a delayed port's flow wraps from the horizon's end — PyPSA's `cyclic_delay`,
      `cyclic_delay2`, …; where it does not, the flow still in transit at the first snapshots is lost
    dims: [link_output]
    dtype: bool
  Link_marginal_cost:
    description: cost of one unit of flow
    dims: [snapshot, link]
  Load_p_set:
    description: demand
    dims: [snapshot, load]
  snapshot_weightings_generators:
    description: PyPSA's `snapshot_weightings.generators` — hours a snapshot stands for in an energy total
    dims: [snapshot]
  Generator_p_nom_min:
    description: least 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]
  Link_p_nom_min:
    description: least nominal power an extendable link may be built at
    dims: [link]
  Link_capital_cost:
    description: cost of one unit of nominal power — PyPSA's `capital_cost`, periodized as an annuity
      in data prep
    dims: [link]
  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_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]
  GlobalConstraint_type:
    description: which formula the row takes — `primary_energy`, `operational_limit`, `transmission_volume_expansion_limit`,
      `transmission_expansion_cost_limit` or `tech_capacity_expansion_limit`
    dims: [global_constraint]
    dtype: str
  GlobalConstraint_sense:
    description: which way the row binds — `<=`, `>=` or `==`
    dims: [global_constraint]
    dtype: str
  GlobalConstraint_constant:
    description: the constant the total is held against; what a variable cannot carry — an initial charge,
      a non-extendable build — is folded in here by data prep
    dims: [global_constraint]
  Generator_primary_energy_weight:
    description: the constrained attribute per unit of energy at the bus — the carrier's `co2_emissions`
      over the generator's efficiency, data prep; a generator of an unweighted carrier has no row
    dims: [global_constraint, generator]
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]
  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`, lossless'
    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
  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
  Link_p_nom_ext:
    description: '`Link-p_nom` — nominal power where it is a decision; the parameter of the same PyPSA
      name carries the fixed regime'
    dims: [link]
    where: Link_p_nom_extendable
constraints:
  Generator_ext_p_lower:
    description: '`Generator-ext-p-lower` — an extendable generator outputs at least its minimum of the
      chosen build'
    dims: [snapshot, generator]
    where: Generator_p_nom_extendable AND not Generator_committable
    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: [snapshot, generator]
    where: Generator_p_nom_extendable AND not Generator_committable
    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
  Link_ext_p_lower:
    description: '`Link-ext-p-lower` — an extendable link carries at least its minimum of the chosen build,
      negative for the other way'
    dims: [snapshot, link]
    where: Link_p_nom_extendable
    expression: Link_p >= Link_p_min_pu * Link_p_nom_ext
  Link_ext_p_upper:
    description: '`Link-ext-p-upper` — an extendable link carries at most the chosen build'
    dims: [snapshot, link]
    where: Link_p_nom_extendable
    expression: Link_p <= Link_p_max_pu * Link_p_nom_ext
  Link_ext_p_nom_lower:
    description: '`Link-ext-p_nom-lower` — the chosen build is at least its floor'
    dims: [link]
    where: Link_p_nom_extendable
    expression: Link_p_nom_ext >= Link_p_nom_min
  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_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_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
  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
  GlobalConstraint_primary_energy_ub:
    description: '`primary_energy` — its total, at most its constant'
    dims: [global_constraint]
    where: GlobalConstraint_type == 'primary_energy' AND GlobalConstraint_sense == '<='
    expression: primary_energy <= GlobalConstraint_constant
  Bus_nodal_balance:
    description: '`Bus-nodal_balance` — what is generated at a bus, storage dispatch and stores included,
      less what the links take away, plus what arrives over them after losses and any delay at every port
      they deliver to, meets the load there. A bus nothing is attached to has no row; PyPSA refuses one
      that carries load, and this file does not yet.'
    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(Link_output_arrival, 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)
      == sum(Load_p_set, by=Load_bus, over=load, into=bus)
expressions:
  Link_output_arrival:
    description: what a link delivers to an output port at a snapshot — its flow after the port's efficiency,
      delayed by the port's `delay`; where the port is `cyclic_delay` the delayed flow wraps from the
      horizon's end, and where it is not the flow still in transit at the first snapshots is lost. A port
      that does not delay (`delay` zero) delivers its flow unshifted, cyclic or not
    dims: [snapshot, link_output]
    cases:
      wrapping: {when: Link_output_cyclic_delay, expression: 'shift(at(Link_p, by=Link_output_link, over=link,
          into=link_output) * Link_efficiency, along=snapshot, offset=Link_output_delay, edge=''wrap'')'}
    otherwise: shift(at(Link_p, by=Link_output_link, over=link, into=link_output) * Link_efficiency, along=snapshot,
      offset=Link_output_delay, edge=0)
  primary_energy: {description: 'what a `primary_energy` row totals — weighted generator energy, less
      the charge left in weighted storage at the horizon''s end; the initial charge it is compared against
      is folded into the row''s constant', expression: 'sum(sum(Generator_p * snapshot_weightings_generators
      * Generator_primary_energy_weight, over=snapshot), over=generator)'}
objective: {sense: minimize, description: 'operating cost, each snapshot weighted by the hours it stands
    for', expression: sum(Generator_p * Generator_marginal_cost * snapshot_weightings_objective) + sum(Link_p
    * Link_marginal_cost * snapshot_weightings_objective) + sum(Generator_p_nom_ext * Generator_capital_cost)
    + sum(Link_p_nom_ext * Link_capital_cost) + 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 _emissions(n: pypsa.Network, gc: pd.Series) -> pd.Series:
    """The nonzero values of the carrier attribute a `primary_energy` row weighs."""
    values = n.carriers[gc['carrier_attribute']]
    return values[values != 0]


def _gc_constants(n: pypsa.Network) -> pd.DataFrame:
    """Each row's constant, net of the initial charge PyPSA folds into its side of the row.

    A `primary_energy` or `operational_limit` row counts what its non-cyclic
    storage draws down, so PyPSA adds the initial charge as a constant on the
    variable side; the file keeps the variables and moves it here.
    """
    rows = []
    for label, gc in n.global_constraints.iterrows():
        constant = float(gc['constant'])
        if gc['type'] == 'primary_energy':
            emissions = _emissions(n, gc)
            sus = n.storage_units
            member = sus['carrier'].isin(emissions.index) & ~sus['cyclic_state_of_charge']
            constant -= float(
                (sus.loc[member, 'carrier'].map(emissions) * sus.loc[member, 'state_of_charge_initial']).sum()
            )
            stores = n.stores
            member = stores['carrier'].isin(emissions.index) & ~stores['e_cyclic']
            constant -= float((stores.loc[member, 'carrier'].map(emissions) * stores.loc[member, 'e_initial']).sum())
        if gc['type'] == 'operational_limit':
            sus = n.storage_units
            member = (sus['carrier'] == gc['carrier_attribute']) & ~sus['cyclic_state_of_charge']
            constant -= float(sus.loc[member, 'state_of_charge_initial'].sum())
            stores = n.stores
            member = (stores['carrier'] == gc['carrier_attribute']) & ~stores['e_cyclic']
            constant -= float(stores.loc[member, 'e_initial'].sum())
        rows.append({'global_constraint': str(label), 'value': constant})
    return pd.DataFrame(rows, columns=['global_constraint', '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),
    'load': pl.Series('load', list(names(loads.index).astype(str)), 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),
    'global_constraint': pl.Series(
            'global_constraint', list(names(n.global_constraints.index).astype(str)), dtype=pl.String
        ),
        **scenarios(n),
        **periods(n),
        **carriers(n),
    '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'),
    'Line_bus0': relation(n, 'Line', 'bus0'),
    'Line_bus1': relation(n, 'Line', 'bus1'),
    'snapshot_weightings_objective': weighting(n, 'objective'),
    'Generator_p_nom_extendable': static(n, 'Generator', 'p_nom_extendable'),
    '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_committable': static(n, 'Generator', 'committable'),
    'Link_p_nom_extendable': static(n, 'Link', 'p_nom_extendable'),
    '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_output_delay': _per_port(n, 'delay', 'value'),
    'Link_output_cyclic_delay': _per_port(n, 'cyclic_delay', 'value'),
    'Link_marginal_cost': varying(n, 'Link', 'marginal_cost'),
    'Load_p_set': varying(n, 'Load', 'p_set'),
    'snapshot_weightings_generators': weighting(n, 'generators'),
    'Generator_p_nom_min': static(n, 'Generator', 'p_nom_min'),
    'Generator_capital_cost': static(n, 'Generator', 'capital_cost'),
    'Link_p_nom_min': static(n, 'Link', 'p_nom_min'),
    'Link_capital_cost': static(n, 'Link', 'capital_cost'),
    '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_capital_cost': static(n, 'Line', 'capital_cost'),
    'Line_cycle_weight': _cycle_weights(n),
    'GlobalConstraint_type': static(n, 'GlobalConstraint', 'type').astype({'value': str}),
    'GlobalConstraint_sense': static(n, 'GlobalConstraint', 'sense').astype({'value': str}),
    'GlobalConstraint_constant': _gc_constants(n),
    'Generator_primary_energy_weight': _weights(
            primary, generators, 'generator', lambda gc, g: _emissions(n, gc).get(g['carrier'], 0.0) / g['efficiency']
        ),
}

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

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

"""Rung 11: PyPSA's own `ac_dc_meshed` example, whole — meshed AC and DC, extendable lines, links and generators, carriers, a CO2 budget."""

from __future__ import annotations

from datetime import datetime

#: Ten hourly stamps, the example's own. Every weighting column there is 1.0, which is
#: also the default, so no row below sets one.
SNAPSHOTS = [datetime(2015, 1, 1, hour) for hour in range(10)]

#: Wind availability per snapshot, for the three generators that carry a profile.
P_MAX_PU = {
    'Manchester Wind': [0.930019875, 0.4857475804, 0.2336917351, 0.2576042221, 0.6269055694, 0.6035984088, 0.6789075462, 0.3613026112, 0.6216040549, 0.5215183715],
    'Norway Wind': [0.9745832033, 0.4812903778, 0.4072258018, 0.5999649628, 0.524468219, 0.0096927054, 0.2204533621, 0.8239185004, 0.5562297265, 0.4394160378],
    'Frankfurt Wind': [0.5590784039, 0.7529103711, 0.1234650887, 0.9666766524, 0.8590078044, 0.5261537924, 0.077893008, 0.0590234716, 0.2485544952, 0.1080601728],
}  # fmt: skip

#: Demand per snapshot, for each of the six loads.
P_SET = {
    'London': [35.7962441027, 976.8245614698, 250.5873120464, 130.7531445827, 151.1001686, 931.857051942, 289.8482871447, 864.3433217147, 689.5772637703, 627.8789859434],
    'Frankfurt': [398.0478469638, 432.4361062425, 379.8039282662, 868.3617642835, 548.7707546221, 828.6652426012, 449.2907519075, 699.1637663734, 915.8667802518, 414.8876464034],
    'Norway': [820.035835936, 854.8340468618, 42.550744351, 647.5482327851, 884.0738733306, 509.0624485516, 595.6079648147, 291.6424496984, 2.1534925491, 760.7401765038],
    'Norwich': [415.4625642653, 262.6061464526, 418.4763531902, 552.9595393098, 218.159858091, 791.9762655836, 531.8706808219, 23.5134667186, 970.0590684572, 0.9248336907],
    'Bremen': [640.0863775411, 703.554333706, 440.8361303183, 612.5763056818, 803.4367808051, 605.4006873582, 641.0905902397, 408.0085411725, 912.2477761646, 898.0530916423],
    'Manchester': [857.5514402011, 750.5996237166, 156.5648760141, 527.8708221189, 83.8977589634, 676.6233193474, 731.1371004827, 553.3448891847, 298.338082262, 768.2905859888],
}  # fmt: skip


def build():
    """The example network, stated as the calls that build it.

    A rung states its data inline, so that the PyPSA model under review is the
    script — ``reference.py`` says so and ``test_pypsa_references.py`` checks
    it. The numbers here are PyPSA's own ``ac_dc_meshed``, which is where this
    rung's published objective comes from; ``reference.py`` pins the version
    they were read at.
    """
    import pypsa

    n = pypsa.Network()
    n.set_snapshots(SNAPSHOTS)
    # Bus
    n.add('Bus', 'London', v_nom=380.0, x=-0.13, y=51.5)
    n.add('Bus', 'Norwich', v_nom=380.0, x=1.3, y=52.6)
    n.add('Bus', 'Norwich DC', v_nom=200.0, x=1.3, y=52.5, carrier='DC')
    n.add('Bus', 'Manchester', v_nom=380.0, x=-2.2, y=53.47)
    n.add('Bus', 'Bremen', v_nom=380.0, x=8.8, y=53.08)
    n.add('Bus', 'Bremen DC', v_nom=200.0, x=8.8, y=52.98, carrier='DC')
    n.add('Bus', 'Frankfurt', v_nom=380.0, x=8.7, y=50.12)
    n.add('Bus', 'Norway', v_nom=380.0, x=10.75, y=60.0)
    n.add('Bus', 'Norway DC', v_nom=200.0, x=10.75, y=60.0, carrier='DC')
    # Carrier
    n.add('Carrier', 'gas', co2_emissions=0.24, color='red')
    n.add('Carrier', 'wind', color='blue')
    n.add('Carrier', 'battery', color='green')
    n.add('Carrier', 'load', color='black')
    n.add('Carrier', 'AC', color='orange')
    n.add('Carrier', 'DC', color='purple')
    # Generator
    n.add(
        'Generator',
        'Manchester Wind',
        bus='Manchester',
        p_nom=80.0,
        p_nom_extendable=True,
        p_nom_min=100.0,
        p_max_pu=P_MAX_PU['Manchester Wind'],
        carrier='wind',
        marginal_cost=0.11,
        capital_cost=2793.6516029328,
    )
    n.add(
        'Generator',
        'Manchester Gas',
        bus='Manchester',
        p_nom=50000.0,
        p_nom_extendable=True,
        carrier='gas',
        marginal_cost=4.5323676307,
        capital_cost=196.6151679691,
        efficiency=0.3500264336,
    )
    n.add(
        'Generator',
        'Norway Wind',
        bus='Norway',
        p_nom=100.0,
        p_nom_extendable=True,
        p_nom_min=100.0,
        p_max_pu=P_MAX_PU['Norway Wind'],
        carrier='wind',
        marginal_cost=0.09,
        capital_cost=2184.3747960912,
    )
    n.add(
        'Generator',
        'Norway Gas',
        bus='Norway',
        p_nom=20000.0,
        p_nom_extendable=True,
        carrier='gas',
        marginal_cost=5.8928445406,
        capital_cost=158.2512497168,
        efficiency=0.3568363832,
    )
    n.add(
        'Generator',
        'Frankfurt Wind',
        bus='Frankfurt',
        p_nom=110.0,
        p_nom_extendable=True,
        p_nom_min=100.0,
        p_max_pu=P_MAX_PU['Frankfurt Wind'],
        carrier='wind',
        marginal_cost=0.1,
        capital_cost=2129.4561224763,
    )
    n.add(
        'Generator',
        'Frankfurt Gas',
        bus='Frankfurt',
        p_nom=80000.0,
        p_nom_extendable=True,
        carrier='gas',
        marginal_cost=4.0863219899,
        capital_cost=102.6769530076,
        efficiency=0.3516658529,
    )
    # Line
    n.add(
        'Line',
        '0',
        bus0='London',
        bus1='Manchester',
        x=0.7968782824,
        s_nom=40000.0,
        s_nom_extendable=True,
        capital_cost=0.1367157553,
        carrier='AC',
    )
    n.add(
        'Line',
        '1',
        bus0='Manchester',
        bus1='Norwich',
        x=0.3915599178,
        s_nom=40000.0,
        s_nom_extendable=True,
        capital_cost=0.1334916779,
        carrier='AC',
    )
    n.add(
        'Line',
        '2',
        bus0='Bremen DC',
        bus1='Norwich DC',
        r=0.2126041927,
        s_nom=40000.0,
        s_nom_extendable=True,
        capital_cost=0.0086734246,
        carrier='AC',
    )
    n.add(
        'Line',
        '3',
        bus0='Norwich DC',
        bus1='Norway DC',
        r=0.4861637504,
        s_nom=40000.0,
        s_nom_extendable=True,
        capital_cost=0.1291260515,
        carrier='AC',
    )
    n.add(
        'Line',
        '4',
        bus0='Norway DC',
        bus1='Bremen DC',
        r=0.4287266497,
        s_nom=40000.0,
        s_nom_extendable=True,
        capital_cost=0.0624298729,
        carrier='AC',
    )
    n.add(
        'Line',
        '5',
        bus0='Norwich',
        bus1='London',
        x=0.2388003463,
        s_nom=40000.0,
        s_nom_extendable=True,
        capital_cost=0.0218524519,
        carrier='AC',
    )
    n.add(
        'Line',
        '6',
        bus0='Bremen',
        bus1='Frankfurt',
        x=0.4,
        s_nom=40000.0,
        s_nom_extendable=True,
        capital_cost=0.2,
        carrier='AC',
    )
    # Link
    n.add(
        'Link',
        'Norwich Converter',
        bus0='Norwich',
        bus1='Norwich DC',
        carrier='DC',
        p_nom=1000.0,
        p_nom_extendable=True,
        p_min_pu=-0.9,
        p_max_pu=0.9,
        capital_cost=0.21,
    )
    n.add(
        'Link',
        'Norway Converter',
        bus0='Norway',
        bus1='Norway DC',
        carrier='DC',
        p_nom=1000.0,
        p_nom_extendable=True,
        p_min_pu=-0.9,
        p_max_pu=0.9,
        capital_cost=0.2,
    )
    n.add(
        'Link',
        'Bremen Converter',
        bus0='Bremen',
        bus1='Bremen DC',
        carrier='DC',
        p_nom=1000.0,
        p_nom_extendable=True,
        p_min_pu=-0.9,
        p_max_pu=0.9,
        capital_cost=0.19,
    )
    n.add(
        'Link',
        'DC link',
        bus0='London',
        bus1='Bremen',
        carrier='DC',
        p_nom=1000.0,
        p_nom_extendable=True,
        p_min_pu=-0.9,
        p_max_pu=0.9,
        capital_cost=0.8765342,
    )
    # Load
    n.add('Load', 'London', bus='London', carrier='load', p_set=P_SET['London'])
    n.add('Load', 'Frankfurt', bus='Frankfurt', carrier='load', p_set=P_SET['Frankfurt'])
    n.add('Load', 'Norway', bus='Norway', carrier='load', p_set=P_SET['Norway'])
    n.add('Load', 'Norwich', bus='Norwich', carrier='load', p_set=P_SET['Norwich'])
    n.add('Load', 'Bremen', bus='Bremen', carrier='load', p_set=P_SET['Bremen'])
    n.add('Load', 'Manchester', bus='Manchester', carrier='load', p_set=P_SET['Manchester'])
    # GlobalConstraint
    n.add('GlobalConstraint', 'co2_limit', sense='<=', constant=1000.0)
    return n

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