StorageUnit as Link and Store Components¶
This example demonstrates how the StorageUnit component can be replaced by the more fundamental Store and Link components, and how their parameters map to each other.
In [2]:
Copied!
import pandas as pd
from numpy.testing import assert_almost_equal, assert_array_almost_equal
import pypsa
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.
In [3]:
Copied!
def replace_storage_unit(n: pypsa.Network, name: str) -> tuple:
"""
Replace the storage unit with `name` 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 optimisation to implement this constraint.
"""
su = n.storage_units.loc[name]
bus_name = f"{su['bus']} {su['carrier']}"
link_1_name = f"{name} converter {su['carrier']} to AC"
link_2_name = f"{name} converter AC to {su['carrier']}"
store_name = f"{name} store {su['carrier']}"
gen_name = f"{name} inflow"
n.add("Bus", bus_name, carrier=su["carrier"])
# dispatch link
n.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
n.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 (
name in n.storage_units_t.state_of_charge_set.columns
and (~pd.isnull(n.storage_units_t.state_of_charge_set[name])).any()
):
e_max_pu = pd.Series(data=1, index=n.snapshots)
e_min_pu = pd.Series(data=0, index=n.snapshots)
non_null = ~pd.isnull(n.storage_units_t.state_of_charge_set[name])
e_max_pu[non_null] = n.storage_units_t.state_of_charge_set[name][non_null]
e_min_pu[non_null] = n.storage_units_t.state_of_charge_set[name][non_null]
else:
e_max_pu = 1
e_min_pu = 0
n.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"],
)
n.add("Carrier", "rain", co2_emissions=0)
# inflow from a variable generator, which can be curtailed (i.e. spilled)
inflow_max = n.storage_units_t.inflow[name].max()
if inflow_max == 0:
inflow_pu = 0
else:
inflow_pu = n.storage_units_t.inflow[name] / inflow_max
n.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(n: pypsa.Network, sns: pd.Index) -> None:
m = n.model
lhs = (
m["Store-e_nom"].at[store_name]
- m["Link-p_nom"].at[link_1_name] * ratio1
)
m.add_constraints(lhs == 0, name="store_fix_1")
lhs = (
m["Store-e_nom"].at[store_name]
- m["Link-p_nom"].at[link_2_name] * ratio2
)
m.add_constraints(lhs == 0, name="store_fix_2")
else:
extra_functionality = None
n.remove("StorageUnit", name)
return bus_name, link_1_name, link_2_name, store_name, gen_name, extra_functionality
def replace_storage_unit(n: pypsa.Network, name: str) -> tuple:
"""
Replace the storage unit with `name` 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 optimisation to implement this constraint.
"""
su = n.storage_units.loc[name]
bus_name = f"{su['bus']} {su['carrier']}"
link_1_name = f"{name} converter {su['carrier']} to AC"
link_2_name = f"{name} converter AC to {su['carrier']}"
store_name = f"{name} store {su['carrier']}"
gen_name = f"{name} inflow"
n.add("Bus", bus_name, carrier=su["carrier"])
# dispatch link
n.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
n.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 (
name in n.storage_units_t.state_of_charge_set.columns
and (~pd.isnull(n.storage_units_t.state_of_charge_set[name])).any()
):
e_max_pu = pd.Series(data=1, index=n.snapshots)
e_min_pu = pd.Series(data=0, index=n.snapshots)
non_null = ~pd.isnull(n.storage_units_t.state_of_charge_set[name])
e_max_pu[non_null] = n.storage_units_t.state_of_charge_set[name][non_null]
e_min_pu[non_null] = n.storage_units_t.state_of_charge_set[name][non_null]
else:
e_max_pu = 1
e_min_pu = 0
n.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"],
)
n.add("Carrier", "rain", co2_emissions=0)
# inflow from a variable generator, which can be curtailed (i.e. spilled)
inflow_max = n.storage_units_t.inflow[name].max()
if inflow_max == 0:
inflow_pu = 0
else:
inflow_pu = n.storage_units_t.inflow[name] / inflow_max
n.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(n: pypsa.Network, sns: pd.Index) -> None:
m = n.model
lhs = (
m["Store-e_nom"].at[store_name]
- m["Link-p_nom"].at[link_1_name] * ratio1
)
m.add_constraints(lhs == 0, name="store_fix_1")
lhs = (
m["Store-e_nom"].at[store_name]
- m["Link-p_nom"].at[link_2_name] * ratio2
)
m.add_constraints(lhs == 0, name="store_fix_2")
else:
extra_functionality = None
n.remove("StorageUnit", name)
return bus_name, link_1_name, link_2_name, store_name, gen_name, extra_functionality
Now, take an existing example which has already been solved as a reference case:
In [4]:
Copied!
n_r = pypsa.examples.storage_hvdc()
n_r.optimize(log_to_console=False)
n_r = pypsa.examples.storage_hvdc()
n_r.optimize(log_to_console=False)
INFO:pypsa.network.io:Retrieving network data from https://github.com/PyPSA/PyPSA/raw/v1.0.5/examples/networks/storage-hvdc/storage-hvdc.nc.
INFO:pypsa.network.io:Imported network 'Storage-HVDC' 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='name')
WARNING:pypsa.consistency:The following links have carriers which are not defined: Index(['TL 0', 'TL 1'], dtype='object', name='name')
WARNING:pypsa.consistency:The following lines have carriers which are not defined: Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='name')
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='name')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.09s
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.
Out[4]:
('ok', 'optimal')
We load it again, but now we replace the StorageUnit components with Store and Link components.
In [5]:
Copied!
n = pypsa.examples.storage_hvdc()
name = "Storage 0"
(
bus_name,
link_1_name,
link_2_name,
store_name,
gen_name,
extra_functionality,
) = replace_storage_unit(n, name)
n.optimize(extra_functionality=extra_functionality)
n = pypsa.examples.storage_hvdc()
name = "Storage 0"
(
bus_name,
link_1_name,
link_2_name,
store_name,
gen_name,
extra_functionality,
) = replace_storage_unit(n, name)
n.optimize(extra_functionality=extra_functionality)
INFO:pypsa.network.io:Imported network 'Storage-HVDC' 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='name')
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='name')
WARNING:pypsa.consistency:The following lines have carriers which are not defined: Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='name')
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='name')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io: Writing time: 0.1s
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
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.
Running HiGHS 1.11.0 (git hash: 364c83a): Copyright (c) 2025 HiGHS under MIT licence terms
LP linopy-problem-8189zmxe has 1144 rows; 504 cols; 2413 nonzeros
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
Dependent equations search running on 82 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
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 name : linopy-problem-8189zmxe
Model status : Optimal
Simplex iterations: 423
Objective value : 1.4670508825e+07
P-D objective error : 7.6179160843e-16
HiGHS run time : 0.01
Writing the solution to /tmp/linopy-solve-06e0cyes.sol
Out[5]:
('ok', 'optimal')
Let's ensure that the conversion was accurate by comparing the results of the original and modified networks.
In [6]:
Copied!
assert_almost_equal(n_r.objective, n.objective, decimal=2)
assert_array_almost_equal(
n_r.storage_units_t.state_of_charge[name],
n.stores_t.e[store_name],
)
assert_array_almost_equal(
n_r.storage_units_t.p[name],
-n.links_t.p1[link_1_name] - n.links_t.p0[link_2_name],
)
assert_array_almost_equal(
n_r.storage_units.at[name, "p_nom_opt"],
n.links.at[link_2_name, "p_nom_opt"],
)
assert_array_almost_equal(
n_r.storage_units.at[name, "p_nom_opt"],
n.links.at[link_1_name, "p_nom_opt"]
* n_r.storage_units.at[name, "efficiency_dispatch"],
)
assert_almost_equal(n_r.objective, n.objective, decimal=2)
assert_array_almost_equal(
n_r.storage_units_t.state_of_charge[name],
n.stores_t.e[store_name],
)
assert_array_almost_equal(
n_r.storage_units_t.p[name],
-n.links_t.p1[link_1_name] - n.links_t.p0[link_2_name],
)
assert_array_almost_equal(
n_r.storage_units.at[name, "p_nom_opt"],
n.links.at[link_2_name, "p_nom_opt"],
)
assert_array_almost_equal(
n_r.storage_units.at[name, "p_nom_opt"],
n.links.at[link_1_name, "p_nom_opt"]
* n_r.storage_units.at[name, "efficiency_dispatch"],
)