Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Replace StorageUnits with fundamental Links and Stores#

This notebook demonstrates how storage units can be replaced by more fundamental components, and how their parameters map to each other.

[1]:

import pandas as pd
from numpy.testing import assert_almost_equal, assert_array_almost_equal

import pypsa

We define two functions we use in the following.

[2]:

def replace_su(network, su_to_replace):
    """
    Replace the storage unit su_to_replace with a bus for the energy
    carrier, two links for the conversion of the energy carrier to and from electricity,
    a store to keep track of the depletion of the energy carrier and its
    CO2 emissions, and a variable generator for the storage inflow.

    Because the energy size and power size are linked in the storage unit by the max_hours,
    extra functionality must be added to the LOPF to implement this constraint.
    """

    su = network.storage_units.loc[su_to_replace]

    bus_name = "{} {}".format(su["bus"], su["carrier"])
    link_1_name = "{} converter {} to AC".format(su_to_replace, su["carrier"])
    link_2_name = "{} converter AC to {}".format(su_to_replace, su["carrier"])
    store_name = "{} store {}".format(su_to_replace, su["carrier"])
    gen_name = f"{su_to_replace} inflow"

    network.add("Bus", bus_name, carrier=su["carrier"])

    # dispatch link
    network.add(
        "Link",
        link_1_name,
        bus0=bus_name,
        bus1=su["bus"],
        capital_cost=su["capital_cost"] * su["efficiency_dispatch"],
        p_nom=su["p_nom"] / su["efficiency_dispatch"],
        p_nom_extendable=su["p_nom_extendable"],
        p_nom_max=su["p_nom_max"] / su["efficiency_dispatch"],
        p_nom_min=su["p_nom_min"] / su["efficiency_dispatch"],
        p_max_pu=su["p_max_pu"],
        marginal_cost=su["marginal_cost"] * su["efficiency_dispatch"],
        efficiency=su["efficiency_dispatch"],
    )

    # store link
    network.add(
        "Link",
        link_2_name,
        bus0=su["bus"],
        bus1=bus_name,
        p_nom=su["p_nom"],
        p_nom_extendable=su["p_nom_extendable"],
        p_nom_max=su["p_nom_max"],
        p_nom_min=su["p_nom_min"],
        p_max_pu=-su["p_min_pu"],
        efficiency=su["efficiency_store"],
    )

    if (
        su_to_replace in network.storage_units_t.state_of_charge_set.columns
        and (
            ~pd.isnull(network.storage_units_t.state_of_charge_set[su_to_replace])
        ).any()
    ):
        e_max_pu = pd.Series(data=1.0, index=network.snapshots)
        e_min_pu = pd.Series(data=0.0, index=network.snapshots)
        non_null = ~pd.isnull(
            network.storage_units_t.state_of_charge_set[su_to_replace]
        )
        e_max_pu[non_null] = network.storage_units_t.state_of_charge_set[su_to_replace][
            non_null
        ]
        e_min_pu[non_null] = network.storage_units_t.state_of_charge_set[su_to_replace][
            non_null
        ]
    else:
        e_max_pu = 1.0
        e_min_pu = 0.0

    network.add(
        "Store",
        store_name,
        bus=bus_name,
        e_nom=su["p_nom"] * su["max_hours"],
        e_nom_min=su["p_nom_min"] / su["efficiency_dispatch"] * su["max_hours"],
        e_nom_max=su["p_nom_max"] / su["efficiency_dispatch"] * su["max_hours"],
        e_nom_extendable=su["p_nom_extendable"],
        e_max_pu=e_max_pu,
        e_min_pu=e_min_pu,
        standing_loss=su["standing_loss"],
        e_cyclic=su["cyclic_state_of_charge"],
        e_initial=su["state_of_charge_initial"],
    )

    network.add("Carrier", "rain", co2_emissions=0.0)

    # inflow from a variable generator, which can be curtailed (i.e. spilled)
    inflow_max = network.storage_units_t.inflow[su_to_replace].max()

    if inflow_max == 0.0:
        inflow_pu = 0.0
    else:
        inflow_pu = network.storage_units_t.inflow[su_to_replace] / inflow_max

    network.add(
        "Generator",
        gen_name,
        bus=bus_name,
        carrier="rain",
        p_nom=inflow_max,
        p_max_pu=inflow_pu,
    )

    if su["p_nom_extendable"]:
        ratio2 = su["max_hours"]
        ratio1 = ratio2 * su["efficiency_dispatch"]

        def extra_functionality(network, snapshots):
            model = network.model
            model.add_constraints(
                model["Store-e_nom"][store_name]
                - model["Link-p_nom"][link_1_name] * ratio1
                == 0,
                name="store_fix_1",
            )
            model.add_constraints(
                model["Store-e_nom"][store_name]
                - model["Link-p_nom"][link_2_name] * ratio2
                == 0,
                name="store_fix_2",
            )

    else:
        extra_functionality = None

    network.remove("StorageUnit", su_to_replace)

    return bus_name, link_1_name, link_2_name, store_name, gen_name, extra_functionality

Now, take an example from the git repo which has already been solved

[3]:

network_r = pypsa.examples.storage_hvdc(from_master=True)
network_r.optimize()

WARNING:pypsa.io:Importing network from PyPSA version v0.17.1 while current version is v0.29.0. Read the release notes at https://pypsa.readthedocs.io/en/latest/release_notes.html to prepare your network for import.
INFO:pypsa.io:Imported network storage-hvdc.nc has buses, carriers, generators, global_constraints, lines, links, loads, storage_units
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['TL 0', 'TL 1'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following lines have carriers which are not defined:
Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='Line')
WARNING:pypsa.consistency:The following lines have zero r, which could break the linear load flow:
Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='Line')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['TL 0', 'TL 1'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following lines have carriers which are not defined:
Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='Line')
WARNING:pypsa.consistency:The following lines have zero r, which could break the linear load flow:
Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='Line')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io: Writing time: 0.18s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 490 primals, 1104 duals
Objective: 1.47e+07
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Line-fix-s-lower, Line-fix-s-upper, Line-ext-s-lower, Line-ext-s-upper, Link-fix-p-lower, Link-fix-p-upper, Link-ext-p-lower, Link-ext-p-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-state_of_charge_set, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [8e-05, 6e+00]
  Cost   [3e-03, 3e+03]
  Bound  [2e+01, 1e+06]
  RHS    [2e+00, 1e+03]
Presolving model
673 rows, 487 cols, 1922 nonzeros  0s
585 rows, 401 cols, 2258 nonzeros  0s
Presolve : Reductions: rows 585(-519); columns 401(-89); elements 2258(-101)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0    -1.4501656367e+05 Ph1: 385(1.3432e+06); Du: 125(265.685) 0s
        413     1.4670508825e+07 Pr: 0(0); Du: 0(9.89209e-14) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 413
Objective value     :  1.4670508825e+07
HiGHS run time      :          0.01
Writing the solution to /tmp/linopy-solve-lp5elq72.sol

[3]:

('ok', 'optimal')

[4]:

network = pypsa.examples.storage_hvdc(from_master=True)

su_to_replace = "Storage 0"

(
    bus_name,
    link_1_name,
    link_2_name,
    store_name,
    gen_name,
    extra_functionality,
) = replace_su(network, su_to_replace)
network.optimize(extra_functionality=extra_functionality)

assert_almost_equal(network_r.objective, network.objective, decimal=2)
assert_array_almost_equal(
    network_r.storage_units_t.state_of_charge[su_to_replace],
    network.stores_t.e[store_name],
)
assert_array_almost_equal(
    network_r.storage_units_t.p[su_to_replace],
    -network.links_t.p1[link_1_name] - network.links_t.p0[link_2_name],
)

# check optimised size
assert_array_almost_equal(
    network_r.storage_units.at[su_to_replace, "p_nom_opt"],
    network.links.at[link_2_name, "p_nom_opt"],
)
assert_array_almost_equal(
    network_r.storage_units.at[su_to_replace, "p_nom_opt"],
    network.links.at[link_1_name, "p_nom_opt"]
    * network_r.storage_units.at[su_to_replace, "efficiency_dispatch"],
)

WARNING:pypsa.io:Importing network from PyPSA version v0.17.1 while current version is v0.29.0. Read the release notes at https://pypsa.readthedocs.io/en/latest/release_notes.html to prepare your network for import.
INFO:pypsa.io:Imported network storage-hvdc.nc has buses, carriers, generators, global_constraints, lines, links, loads, storage_units
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['TL 0', 'TL 1', 'Storage 0 converter AC to battery'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following lines have carriers which are not defined:
Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='Line')
WARNING:pypsa.consistency:The following lines have zero r, which could break the linear load flow:
Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='Line')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['TL 0', 'TL 1', 'Storage 0 converter AC to battery'], dtype='object', name='Link')
WARNING:pypsa.consistency:The following lines have carriers which are not defined:
Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='Line')
WARNING:pypsa.consistency:The following lines have zero r, which could break the linear load flow:
Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='Line')
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.29.0/lib/python3.12/site-packages/linopy/variables.py:198: FutureWarning: Accessing a single value with `Variable[...]` and return type ScalarVariable is deprecated. In future, this will return a Variable.To get a ScalarVariable use `Variable.at[...]` instead.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io: Writing time: 0.21s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 504 primals, 1144 duals
Objective: 1.47e+07
Solver model: available
Solver message: optimal

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [8e-05, 6e+00]
  Cost   [3e-03, 3e+03]
  Bound  [1e+06, 1e+06]
  RHS    [2e+00, 1e+03]
Presolving model
682 rows, 498 cols, 1942 nonzeros  0s
585 rows, 401 cols, 2254 nonzeros  0s
585 rows, 401 cols, 2254 nonzeros  0s
Presolve : Reductions: rows 585(-559); columns 401(-103); elements 2254(-159)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0    -1.4523089525e+05 Ph1: 390(1.38972e+06); Du: 128(270.398) 0s
        423     1.4670508825e+07 Pr: 0(0); Du: 0(7.22311e-13) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 423
Objective value     :  1.4670508825e+07
HiGHS run time      :          0.01
Writing the solution to /tmp/linopy-solve-r1fdcwnu.sol

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Line-fix-s-lower, Line-fix-s-upper, Line-ext-s-lower, Line-ext-s-upper, Link-fix-p-lower, Link-fix-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-state_of_charge_set, Kirchhoff-Voltage-Law, StorageUnit-energy_balance, Store-energy_balance, store_fix_1, store_fix_2 were not assigned to the network.